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Abstract 

We provide the explicit solutions of linear, left-invariant, (convection)-diffusion equations 
and the corresponding resolvent equations on the 2D-Euclidean motion group SE(2) — R 2 xi T. 
These diffusion equations are forward Kolmogorov equations for well-known stochastic pro- 
cesses for contour enhancement and contour completion. The solutions are given by group- 
convolution with the corresponding Green's functions which we derive in explicit form. We 
have solved the Kolmogorov equations for stochastic processes on contour completion, in 
earlier work [T5] . Here we mainly focus on the Forward Kolmogorov equations for contour en- 
hancement processes which, in contrast to the Kolmogorov equations for contour completion, 
do not include convection. The Green's functions of these left-invariant partial differential 
equations coincide with the heat-kernels on SE(2). Nevertheless, our exact formulae do not 
seem to appear in literature. Furthermore, by approximating the left-invariant basis of the 
generators on SE(2) by left-invariant generators of a Heisenberg-type group, we derive ap- 
proximations of the Green's functions. 

The Green's functions are used in so-called completion distributions on SE(2) which are 
the product of a forward resolvent evolved from a source distribution on SE(2) and a backward 
resolvent evolution evolved from a sink distribution on SE(2). Such completion distributions 
on SE(2) represent the probability density that a random walker from a forward proces collides 
with a random walker from a backward process. On the one hand, the modes of Mumford's 
direction process (for contour completion) coincides with elastica curves minimizing J K 2 + eds, 
and they are closely related to zero- crossings of two left-invariant derivatives of the comple- 
tion distribution. On the other hand, the completion measure for the contour enhancement 
proposed by Citti and Sarti, concentrates on the geodesies minimizing J \J n? + eds if the 
expected life time 1/a of a random walker in SE(2) tends to zero. 

This motivates a comparison between the geodesies and elastica. For reasonable parameter 
settings they turn out to be quite similar. However, we apply the results by Bryant and 
Griffiths 9 on Marsden-Weinstein reduction on Euler-Lagrange equations associated to the 
elastica functional, to the case of the geodesic functional. This yields rather simple practical 
analytic solutions for the geodesies, which in contrast to the formula for the elastica, do not 
involve special functions. 

The theory is directly motivated by several medical image analysis applications where 
enhancement of elongated structures, such as catheters and bloodvessels, in noisy medical 
image data is required. Within this article we show how the left-invariant evolution processes 



can be used for automated contour enhancement/completion using a so-called orientation 
score, which is obtained from a grey-value image by means of a special type of unitary wavelet 
transformation. Here the (invertible) orientation score serves as both the source and sink- 
distribution in the completion distribution. 

Furthermore, we also consider non-linear adaptive evolution equations on orientation scores. 
These non-linear evolution equations are practical improvements of the celebrated standard 
"coherence enhancing diffusion" -schemes on images as they can cope with crossing contours. 
Here we employ differential geometry on SE(2) to include curvature in our non-linear diffusion 
scheme on orientation scores. Finally, we use the same differential geometry for a morphology 
theory on orientation scores yielding automated erosion towards geodesics/elastica. 
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1 Invertible Orientation Scores 

In many image analysis applications an object Uf <E L 2 (SE(2) ) defined on the 2D-Euclidean 
motion group SE{2) = K 2 xi T is constructed from a 2D-grey- value image / € L 2 (IR 2 ). Such an 
object provides an overview of all local orientations in an image. This is important for image 
analysis and perceptual organization, [35], [25], [U], [23], [2D], [55], [7] and is inspired by our own 
visual system, in which receptive fields exist that are tuned to various locations and orientations, 
[51] , [8j. In addition to the approach given in the introduction other schemes to construct Uf : 
M 2 x: T — > C from an image / : M 2 — > M. exist, but only few methods put emphasis on the stability 
of the inverse transformation Ut 1— * /. 

In this section we provide an example on how to obtain such an object Uf from an image /. 
This leads to the concept of invertible orientation scores, which we developed in previous work, 
[15] . [20] . [18] . and which we briefly explain here. 

An orientation score Uf : R 2 X T — > C of an image / : E 2 — * K is obtained by means of an 
anisotropic convolution kernel ip : R 2 — » C via 

U f (g) = [ V(^ 1 (y-x))/(y) dy, g = (x, e* 9 ) e G = M 2 x T, R e G SO(2), 

where ip(—x) — V'(x). Assume ^ € L 2 (R 2 ) n L^M 2 ), then the transform which maps image 
/ G L 2 (IR 2 ) onto its orientation score Uf £ L2(M 2 xi T) can be re- written as 

Uf{g) = {W i ,f){g) = {U g ^f)^ m , 

where g 1— > W ff is a unitary (group-)representation of the Euclidean motion group SE{2) = M 2 xi T 
into L 2 (]R 2 ) given by U g f(y) = /{R^iy - x)) for all 3 = (x,e 4e ) S S , J B(2) and all / e L 2 (E 2 ). 
Note that the representation U is reducible as it leaves the following closed subspaces invariant 
{/ G L 2 (IR 2 ) I supp^f/] C B Q e }, g > 0, where B , e denotes the ball with center Oel 2 and 
radius g > and where T : L 2 (K 2 ) — > L 2 (M 2 ) denotes the Fourier transform given by 

•^/M = ^ | /(x) e -^> x >dx, 

R 2 

for almost every u) e R 2 and all / e L 2 (M 2 ). 

This differs from standard continuous wavelet theory, see for example [39] and [4], where the 
wavelet transform is constructed by means of a quasi-regular representation of the similitude 
group M. d xi T x M + , which is unitary, irreducible and square integrable (admitting the application 
of the more general results in [33]). For the image analysis this means that we do allow a stable 
reconstruction already at a single scale orientation score for a proper choice of ip. In standard 
wavelet reconstruction schemes, however, it is not possible to obtain an image / in a well-posed 
manner from a "fixed scale layer", that is from W^/(-, •,&) 6 L 2 (M 2 xi T), for fixed scale a > 0. 

Moreover, the general wavelet reconstruction results [33J do not apply to the transform / 1—* Uf, 
since our representation hi is reducible. In earlier work we therefore provided a general theory 
[T5] . [T2"] . [13] . to construct wavelet transforms associated with admissible vectors/ distributions^] 
With these wavelet transforms we construct orientation scores Uf : K 2 xi T — ► C by means of 

1 Depending whether images are assumed to be band-limited or not, for full details see |14j . 
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admissible line detecting vectors 2 ] ip <G L2(K 2 ) such that the transform is unitary onto the 
unique reproducing kernel Hilbert space C^- of functions on SE(2) with reproducing kernel 
K(g, h) = (U g tp,Uhip), which is a closed vector subspace of h 2 (SE(2)). For the abstract construc- 
tion of the unique reproducing kernel space <C l K on a set I (not necessarily a group) from a function 
of positive type K : I x I — > C, we refer to the early work of Aronszajn [5]. Here we only provide 
the essential Planchercl formula, which can also be found in a slightly different way in the work 
of Fiihr [3D], for the wavelet transform and which provides a more tangible description of the 
norm on cff^ rather than the abstract one in [S]. To this end we note that we can write 



(nV)(x,e* e ) = (W (x , e ^,/) L2(R2) = {TT^lZ e ^,Tf)^ m = ^(llgFiP ■ ^/)(x) 

where the rotation and translation operators on L 2 (IR 2 ) are defined by 7Zgf(y) — f(R e ~ 1 y) and 
7x/(y) = /(y — x). Consequently, we find that 

imfWl^ =//|(JTV*/)(w,e <<, )| a <W B ^ 6j jdu» 

^})\{Tf)(^\T^^A6 w ^- ) Au> (1.1) 



R 2 T 



J R2 |(^/)(u,)| 2 du = \\f\\l 2{t 



■V 



where € C(K 2 ,K) is given by M^(w) = \Tip(Rj(jj)\ 2 d9. If ip is chosen such that M = 1 
then we gain L,2-norm preservation. However, this is not possible as ip G ]L.2(R 2 ) nLi(R 2 ) implies 
that Mif, is a continuous function vanishing at infinity. Now theoretically speaking one can use 
a Gelfand-triple structure generated by ^/l + |A| to allow distributional wavelets^] ip E H _fe (R 2 ), 
k > 1 with the property = 1, so that ip has equal length in each irreducible subspace (which 
uniquely correspond to the dual orbits of SO(2) on K 2 ), for details and generalizations see [2]. In 
practice, however, because of finite grid sampling, we can as well restrict hi (which is well-defined) 
to the space of bandlimited images. 

Finally, since the wavelet transform maps the space of images L<2(K 2 ) unitarily onto the 

S E( 2) 

space of orientation scores C K (provided that Mj, > 0) we can reconstruct the original image 
/ : R 2 — > K from its orientation score Uf : SE{2) — » C by means of the adjoint 



/ = W;W V 4/] = I 7 ' 1 [w » Jo T[U f (;e^)](u>) F[Ke^)H d8 M^(u>)\ (1.2) 

For typical examples (and different classes) of wavelets ip such that = 1 and details on fast 
approximative reconstructions see [55], [T5],[T7]. For an illustration of a typical proper wavelet 
ip (i.e. M.0 « 1) with corresponding transformation W^f and corresponding M$ : R 2 — ► R + 
usually looks like in our relatively fast algorithms, working with discrete subgroups of the torus, 
see Figure [l] With this well-posed, unitary transformation between the space of images and the 
space of orientation scores at hand, we can perform image processing via orientation scores, see 
[TT] , [T5] . [17] . [2U] . [55] , However, for the remainder of the article we assume that the object Uf 

Q fP(n\ 

is some given function in h 2 (SE(2)) and we write U € h 2 (SE(2)) rather than U f € C K v '. For 
all image analysis applications where an object Ut S h 2 (SE(2)) is constructed from an image 
/ € L 2 (R 2 ), operators on the object U £ h 2 (SE(2)) must be left-invariant to ensure Euclidean 
invariant image processing [15jp. 153. This applies also to the cases where the original image cannot 
be reconstructed in a stable manner as in channel representations |25j and steerable tensor voting 



2 Left-invariant Diffusion on the Euclidean Motion Group 

The group product within the group SE(2) of planar translations and rotations is given by 
gg' = (x,e ie )(x',e* e ') = (x + i?,x', e^+ e ')), g = {x,e ie ),g' = (x',e^') E SE(2), 

2 Or rather admissible distributions ij) S H _ ' 1+e '> 2 (IR 2 ), e > if one does not want a restriction to bandlimited 
images. 

^Just like the Fourier transform on L 2 (IR 2 ), where x >-> e iw,x is not within L 2 (M 2 ). 
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Figure 1: (a) Example image (x,y) i— > f(x,y). (b) The structure of the corresponding ori- 
entation score Uf := W^[/]. The circles become spirals and all spirals are situated in the 
same helicoid-shaped plane, (c) Real part of orientation score (x,y) i— > Uf(x, y,e t6 ) displayed 
for 4 different fixed orientations, (d) The absolute vale (x,y) i— > \Uf(x,y,e )\ yields a phase- 
invariant and positive response displayed for 4 fixed orientations, (e) Real part of the wavelet 

||x|| 2 , o 2 

I n 9 ((</>mod2ir)-f ) 



^(X) = 



, _ . , M(p)](x), where M(o) = 

I and Nyquist frequency g and fc-th order B-spline Bk = Bq * k Bq and Bq(x) 



lr 



with 

and parameter values k — 2, g = 4, |er 2 = 400, s = 10, rig = 64. (f) Imaginary part of ip. (g) The 
function | 2 (h) The function M^. 



with Ra = 



G 5*0(2). The tangent space at the unity element e = (0,0, e M ), 



cos ft — sin ( 
sin 8 cos 

T e (SE(2)), is a 3D Lie algebra equipped with Lie product L4,B] = lim uo (a(t)&(t)(a(t)) _1 (&(t)) -1 - e) 
where t i— > a(i) resp. £ i— > 6(t) are any smooth curves in G with a(0) = 6(0) = e and a'(0) = A 
and 6'(Q) = B. Define A 2 , A 3 } := {e 6l e x , e y }. Then {A U A 2 ,A 5 } form a basis of T e (SE(2)) 
and their Lie-products are 



[A 1 ,A 2 ] = A 3 , [A U A 3 ] = -A 2 , [A 2 ,A 3 } = 



(2.3) 



A vector field on SE(2) is called left-invariant if for all g <E G the push-forward of (L g )^X e by left 
multiplication L g h — gh equals X g , that is 



(X g ) = (L g ),(X e ) * fl / - X e (/ o L 9 ), for all / e C° 



(2.4) 



where Q g is some open set around g G SE(2). Recall that the tangent space at the unity element 
e = (0,0, e i0 ), T e (G), is spanned by T e {G) = span{e e , e*, e y } = span{(l, 0, 0), (0, 1, 0), (0, 0, 1)}. By 
the general recipe of constructing left-invariant vector fields from elements in the Lie-algebra 
T e (G) (via the derivative of the right regular representation) we get the following basis for the 
space £(SE(2)) of left-invariant vector fields : 

{-4i,_4 2 , A-s} = {do,d£,d v } = {dg,cos9d x + sm9d y , — sm.6d x + cos9d y }, (2.5) 

with £ = x cos 9 + y sm9, r\ — —x sm6 + y cos 9. More precisely, the left-invariant vector-fields 
are given by 



e#(x, e l9 ) — eg, e^(x, e 10 ) — cos9 e x +sm9 e y , e v (x.,e" 



-sin 9 e^ + cos 9 e y , 



(2.6) 
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Left invariance of tangent vectors to curves 




Figure 2: Left invariant vector fields on SE(2), where we both consider the tangent vectors tangent 
to curves, that is X g = c 1 eg(g) + c 2 e^(g) + c 3 e v (g) for all g £ SE(2), and as differential operators 
on locally defined smooth functions, that is X g = c 1 dg\ + c 2 + c 3 d n \ for all g £ SE(2). 



We see that the push forward of the left multiplication connects the tangent space T e {SE(2)) 
to all tangent spaces T g (SE(2)) . Conversely, the Cartan connection D = d + u> (6.591 on the 
vector bundle (SE(2) 7 T(SE(2))) connects all tangent spaces to T e (SE(2)). Here we note that 



7 (X g ) := (VO, 



Xg - X ( 



for all left-invariant vector fields X. 



where we identified T g={x ^ e) (R 2 , e ie ) with T e (R 2 ,e l °) and T g=(x ^e } (x,T) with T e (0,T), by par- 
allel transport (on K 2 respectively T). We can always consider these vector fields as differential 



operators (i.e. replace by di, i — 9,£,,rf), which yields (2.5). Summarizing, we see that for 
left-invariant vector fields the tangent vector at g is related to the tangent vector at e by ( |2.4[ ) . In 
fact, the push forward (L g )* of the left multiplication puts a Cartan-connectiorj^] between tangent 



spaces, T e (SE(2)) and T g (SE{2)). Equality §2M sets the isomorphism between T e {SE{2)) and 
C{SE(2)), as At <-> A t , i = 1,2,3 implies [Ai,A~\ <-> [Ai,Aj], j = 1,2,3, recall fj2^3}. Moreover it 
is easily verified that 

[Ai,A 2 ] = A 1 A 2 -A 2 A 1 = A 3 , [A 1 ,A 3 ] = -A 2 , [A 2 ,A 3 } = 0. 

See Figure |2] for a geometric explanation of left invariant vector fields, both considered as tangent 
vectors to curves in SE{2) and as differential operators on locally defined smooth functions. 
Example: 



4 This Cartan connection can be related to a left-invariant metric induced by the Killing-form, which is degen- 
erate on SE(2). This can be resolved by pertubing the vectorfields into C(SO(3)) = so(3) by {— f3 2 y COS 9 dg + 
cos 69 ^l + P 2 y 2 d x +sm9(l+(3 2 y 2 ) d v , f3 2 y sin 6 dg-sin 6 s/l+P 2 y 2 d x +cos 9 (1+/3 2 y 2 ) d v , dg}, < /3 « 1. 



See subsection 6.1 
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Consider e x s T e (G), then the derivative of the right-regular representation gives us 



(d^(A 1 )$)( ff ) = (dU(e x )^)(g) = lim 

hlO 

= lim 



4>(g e he 



-<&(9) 



lim 



$(x+/t(cosfl,sine),e' a )-<t'(x,e' £ ') 



-4i*(ff) 



(2.7) 



(cos (9 a, + sin6»5 y )$(g) = <%<%), 



for all $ smooth and defined on some open environment around g = (x, e 10 ) € G. 

Next we follow our general theory for left-invariant scale spaces on Lie-groups, see [16], and 
set the following quadratic form on £(SE(2)) 



(2.8) 



Q D ' a (^li,^2,^3) = ]T -a,A + ^ DijAiAj , Oi, G K,Z> := [Aj] > 0, D T = D 
i=l V 3=1 / 

and consider the only linear left-invariant 2nd-order evolution equations 



d s W = Q a (A 1 ,A2,A 3 ) W 
limW(;s) = U f (-) . 

sj.0 



with corresponding resolvent equations (obtained by Laplace transform over s): 



P = a(Q°^(A 1 ,A2,A 3 ) - aI)~ l Uf. 



(2.9) 



(2.10) 



These resolvent equations are highly relevant as (for the cases a = 0) they correspond to first 
order Tikhonov regularizations on SE(2), |16j . |llj . They also have an important probabilistic 
interpretation, as we will explain next. 

By the results in [15], [33], [T5], the solutions of these left-invariant evolution equations are 
given by S'i?(2)-convolution with the corresponding Green's function: 



W(g,s) = (G»-** S E(2)U)(g)= J G^' a (h~ 1 g)U(h) dfx S E(2){h), g(x,e id ), 

SE(2) 

= I / G?' a (i? e ; 1 (x-x'),e l ( e - e '))i7 / (x',e i0 ')df'dx' 

K 2 

00 

P a (g) = (R%>» * S£(2) U)(g), R»>» = aj G°.*e— ds. 



(2.11) 



For Gaussian estimates of the Green's functions see the general results in [3S] and [53]. See 
Appendix [D] for details on sharp Gaussian estimates for the Green's functions and formal proof of 
(2.11 1 in the particular case D\\ — D22 > 0, U33 = and a = (which is the Forward Kolmogorov 
equation of the contour enhancement process which we will explain next). 

In the special case Dij = SuSji, a = (k , 1,0) our evolution equation (2.9 1 is the Kolmogorov 
equation 



d s W(g, s) = (^ + D lx d 2 e )W{g, a), g G SE(2), s > 
W(g,0) = U(g) 

of Mumford's direction process, [35] , 

X(s) = X(s) e x + Y(s) e y = X(0) + J Q S cos 6(r) e x + sin 6(r) e y dr, 
6(s) = 6(0) + ^e e + sk ee ~ M(0,2D n ), 



(2.12) 



(2.13) 



for contour completion. The explicit solutions of which we have derived in [19 . 

However, within this article we will mainly focus on stochastic processes for contour enhance- 
ment. For contour enhancement we consider the particular case a = 0. D%s = D31 = D23 = 
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D32 = 0. In particular we consider the case Dij = Sij, D33 = 0, a = so that our evolution 
equation (12. 91), becomes 



d s W(g,s) = (Dnid^ + D^d^W^s) 

W(g,Q) = U(g) ^' i4j 
which is the Kolmogorov equation of the following stochastic process for contour enhancement: 
X(s) = X(0) + y/s e$ J Q S cos 9(r) e x + sin 6(r) e y dr, , . 

e( s ) = e(o) + y?e e , 

with e 4 - J\f(0, 2D 22 ) and e e - AT(0, 2D n ), D n ,D 22 > 0. 

In general the evolution equations ( |2.9[ ) are the forward Kolmogorov equations of all linear left- 
invariant stochastic processes on SE(2), as explained in [TI5], [53]. With respect to this connection 
to probability theory we note that W(g, s) represents the probability density of finding oriented 
random walkei]^] (traveling with unit speed, which allows us to identify traveling time with arc- 
length s) at position g given the initial distribution W(-, 0) — Uf a traveling time s > 0, whereas 
P(g) represents the unconditional probability density of finding an oriented random walker at 
position g given the initial distribution W{-, 0) = Uf regardless its traveling time. To this end we 
note that traveling time T in a Markov process is negatively exponentially distributed 

P(T = s) = ae~ as , 

since this is the only continuous memoryless distribution and indeed a simple calculation yields: 
P(x, y,9\U and T = s) = (Gf « * SE{2) U)(x, y, 9) 

P{x,y,e\ U) = J °° P(x, y, 6 I U and T = s)P{T = s)ds = (R^ * SE{2) U)(x,y,9) (2.16) 
with i?f 11 = a J R+ Gf 11 e" QS ds, 



For exact solutions for the resolvent equations (2.10 1 (in the special case of Mumford's direction 
process), approximations and their relation to fast numerical algorithms, see |19j . 



3 Image Enhancement via left-invariant Evolution Equa- 
tions on Invertible Orientation Scores 

Now that we have constructed a stable transformation between images / and corresponding ori- 
entation scores Uf, in Section [T] we can relate operators T on images to operators $ on orientation 
scores in a robust manner, see Figure [4] It is easily verified that o U g = C g o for all 
g G SE(2), where the left-representation C : G -> B(h 2 (SE(2))) is given by C g $(h) = ^{g^h). 
Consequently, the net operator on the image T is Euclidean invariant if and only if the operator 
on the orientation score is left-invariant, i.e. 



T o U g = U g o T for all g £ SE(2) <S> $ o C g = C g o $ for all g £ SE(2), (3.17) 



see [TSjThm. 21 p. 153. 

Here the diffusions discussed in the previous section, section [2] can be used to construct suitable 
operator $ on the orientation scores. At first glance the diffusions themselves (with certain 
stopping time t > 0) or their resolvents (with parameter a > 0) seem suitable candidates for 
operators on orientation scores, as they follow from stochastic processes for contour enhancement 



5 That is a random walker in the space SE(2) where it is only allowed to move along horizontal curves which 
are curves whose tangent vectors always lie in span{8g,9 f} w hich is the horizontal subspace if we apply the Cartan 
connection on Py = (SE(2), SE(2)/Y, it, R) see section |6. 3 1 In previous work in the field of image analysis, [18] . 
I17| , we called these random walkers "oriented gray value particles" . 
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Figure 3: Top left: six random walks in SE(2) = M 2 xi T (and their projection on M 2 ) of direction 
processes for contour-completion by Mumford (43] with a = (kq, 1, 0), D = diag{Z?n, D22, -D33} for 
various parameter settings of n > and Da > 0. Top middle one random walk (500 steps, with 
step-size 0.005) and its projection to the image plane of the linear left-invariant stochastic process 
for contour enhancement within SE(2) with parameter settings Dn = D22 = \ and -D33 = 
(corresponding to Citti and Sarti's cortical model for contour enhancement, 111]*). Top right; one 
random walk (800 steps, with step-size 0.005) of the stochastic process with parameter settings 

0.75, 



£>n = h&a, D22 — ho' 2 , -D33 = he 2 ,, with ag — 0.75, = 1, a v = 0.5 (other parameters have been 



Bottom: 30 random walks in SE(2) = Mr x T again with Di 



set to zero). 

(7£ = 1, <7 n — 0.5, viewed along 0-axis (left) along x-axis (middle) along y-axis (right). Appropriate 
averaging of infinitely many of these sample paths yields the Green's functions, see Figure [5] of 
the forward Kolmogorov equations (2.9 1. Furthermore we note that Mumford's direction process 
is the only linear left-invariant stochastic process on SE(2) whose sample path projections on 
the image plane are diffcrcntiable. For contour completion this may be a reason to discard the 
other linear left-invariant stochastic processes, see [53]. However, the Green's function of all linear 
left-invariant processes (so also the ones for contour-enhancement) are infinitely diffcrcntiable on 
SE(2) \ {e} iff the Hormander condition as we will discuss in section 4.2.1 see (4.41 ), is satisfied. 
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and contour completion and they even map the space of orientation scores into the space 

of orientation scores cff^ again. But appearances are deceptive since if the operator <!> is left- 
invariant (which must be required, see Figure [4]) and linear then the netto operator T is translation 
and rotation invariant boiling down to an isotropic convolution on the original image, which is of 
course not desirable. 

So our operator <I> must be left-invariant and non-linear and still we would like to directly 
relate such operator to stochastic processes on SE(2) discussed in the previous section. Therefor 
we consider the operators 

V) = ap((Q D ' a (A) - aI)-\UY ((Q D ' & (A))* - al)-\vy)^, (g lg) 

where U (the source distribution) and V (the sink distribution) denote two initial distributions on 
SE(2) and where we take the p-th power of both real part 5ft(f7) and imaginary part 9 ! (J7) sepa- 
rately in a sign-preserving manner, i.e. (U) p means sign{$R(£/)} |!ft([/)| p + i sign{$3(£/)} |S3(E/)| P . 
Here the function $(Z7, V) 6 h 2 (SE(2)) can be considered as the completion distribution^] ob- 
tained from collision of the forwardly evolving source distribution U and backwardly evolving sink 
distribution V, similar to [7]. 

Within this manuscript we shall restrict ourselves to the case where both source and sink equal 
the orientation score of original image /, i.e. U = V = Uf := W^f and only occasionally (for 
example section [5]) we shall study the case where U = S go and V = S gi , where go an d gi are some 
given elements in SE(2). 

In sectionJSJwe shall consider more sophisticated and more practical alternatives to the operator 



given by (3.181. But for the moment we restrict ourselves to the case (3.181 as this is much easier 
to analyse and also much easier to implement as it requires two group convolutions (recall (2.11)) 
with the corresponding Green's functions which we shall explicitly derive in the next section. 

The relation between image and orientation score remains 1-to 1 if we ensure that the operator 
on the orientation score again provides an orientation score of an image: Let C^- B( - 2) denoted the 
space of orientation scores within h2(SE(2)), then the relation is 1-to 1 iff $ maps C^^Mnto 
C^ E ' 2 '. However, we naturally extend the reconstruction to l,2(SE(2)): 

mr'Uig) =T~ X [« ~ f^F[U(;e ie )](u) T[ll^]{") M M~\u)] , (3.19) 

for all U E h2(SE(2)). So the effective part of a operator $ on an orientation score is in fact 



$ where = yV^(}V^) ext is the ortho gonal projection of h 2 (SE(2)) onto C S K E{2) . Recall that 



$ must be left-invariant because of (3.171 



It is not difficult to show that the only linear left-invariant kernel operators on I^2(SE(2)) are 
given by S'i?(2)-convolutions. Recall that these kernel operators are given by ( 2.11| ). Even these 



SE(^) 

S^E(2)-convolutions do not leave the space of orientation scores C K invariant. Although 
{K* SE (2) W^/XflO = / (^h^J)h 2 (M^K(h- 1 g)dn SE (2){h) 

SE{2) 

= ( / U h ipK(h~ 1 g)dii SE (2)(h)J) h2{R 2 ) 

SE(2) 

= ( / Ugh-^ K fr)&VSE(2){h) ,/)l 2 (R2) = (W s ^,/)l 2 (R2) = 
SE{2) 

for all / e L 2 (M 2 ), g £ SE(2), where ip = J se m U^-iip K{h) d[isE(2)(h), the reproducing kernel 

space associated to ip will in general not coincide with the reproducing kernel space associated to ip. 
Here we recall from [13], [15], that ip determines the reproducing kernel K(g, h) = (UgipjUhip)^^)- 

6 In image analysis these distributions are called "completion fields", where the word field is inappropriate. 
7 We use this notation since the space of orientation scores generated by proper wavelet ip is the unique repro- 
ducing kernel space on SE(2) with reproducing kernel K{g,h) = (U g %p Mhi>) , [E]p.221-222, p.120-122 
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Image 

/ e L 2 (R 2 ) 



T V 



Orientation Score 

U f e C S K E{2) C h 2 (SE(2)) 



Processed Image 

t[/] = w; MUf]] 



(w;r xt = 



Processed Score 

$[U f ] e U(SE(2)) 



Contour enhancement 



input image / output image T/ = W1WV$} 




Contour completion 

input image / Output image T/ = MfWfy/ 



-*.-_« • -v. / . 

* f * * * i - - ~ - K S // i K> " 



Figure 4: Top Row: The complete scheme; for admissible vectors if) the linear map W^, is unitary 
from L2(M 2 ) onto a closed subspace C^- B ' 2 ' of ~L,2(SE(2)). So we can uniquely relate a transformation 
on an orientation score to a transformation on an image T = (W^,) ea: o$o W,/>, where 



$ : C 



SE(2) 
K 



c 



SE[2) 
K 



(W^) ext is given by ( 3.19 1. Here we take $ as a concatenation of non-linear invertible greyvalue transforms 



and linear left-invariant evolutions ( 3.18 1, although we stress that for most practical applications it is better 
to replace the operator $ by the adaptive evolution operator W^f ► u(x,y,e 10 ,t) defined by the non- 



linear adaptive left-invariant evolution equation (8.1131 with certain stopping time t > 0. Bottom row 



automated contour enhancement (left) and completion (right). Parameter settings, left; Dij = DaSi 



0, a = 0, right D i} = <5ii<Sji, a = (0, 1,0), 



Dll 



0.1. 
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4 The Heat-Kernels on SE(2). 



In section 4.1 we present the exact formulae, which do not seem to appear in literature, of the 
Green's functions and their resolvents for linear anisotropic diffusion on the group SE(2). Al- 
though the exact resolvent diffusion kernels (which take care of Tikhonov regularization on SE(2), 



[HJ| ) are expressed in only 4 Mathieu functions, we also derive, in section 4.2 the corresponding 
Hcisenberg approximation resolvent diffusion kernels (which are rather Green's functions on the 
space of positions and velocities rather than Green's functions on the space of positions and orien- 
tations) which arise by replacing cos# by 1 and sin# by 9. Although these approximation Green's 
functions are not as simple as in the contour-completion case, [19]ch:4.3, they are more suitable 
if it comes to fast implementations, in particular for the Green's functions of the time processes. 
For comparison between the exact resolvent heat kernels and their approximations, see figure [5] 

4.1 The Exact Heat-Kernels on SE{2) = M 2 x 50(2). 

In this section we will derive the heat-kernels : SE(2) — * IR + and the corresponding resolvent 
kernels R a .o ■ SE(2) — > R + on SE(2). Recall that Si£(2)-convolutk>n with these kernels, see 



(2.11l, provide the solutions of the Forward Kolmogorov equations (2.141 and recall that R a r> 



a J K®e as ds. During this chapter we set D as a constant diagonal matrix. Although U33 = 




(as in (2.14l) has our main interest we also consider the more general case where -D33 > 0. 

The kernels Kf and R a . d are the unique solutions of the respectively the following problems 

f {-D 11 (do) 2 -D 22 (d i ) 2 ~D 33 (d v ) 2 + a)R™ D =aS e , f d,K? = {D lx {d e ) 2 + D 22 (d^) 2 + D 33 (d v ) 2 ) K? 

{ i*£ D (v,0)=.R£ D (v,27r) J VmK?=6 e 

{ RZn G U(SE(2)). { k d £ U{SE{2)) 

The first step here is to perform a Fourier transform with respect to the spatial part = M 2 of 
SE{2) = M 2 xi T , so that we obtain R a<D ,K° G h 2 (SE(2)) n C{SE{2)) given by 

Kf ( Wl , uj 2 ,d) = T[K f (•, -,6)] (ui , wa) . 

R a , D {u)l,U) 2 ,Q) = F[R a ,D{-, ;0)]{u>x,ljJ2). 

Then R a ,D an d satisfy 

(al - Bu)R a ,D = %6o nndd s K? =B U K°, limiff (w, 0) = <J e (4.20) 

slO v ' 

where we define the operator 

= -D 22P 2 cos 2 (<^ -6)- D 33 p 2 sin 2 (cp -9) + D n (d e ) 2 

where we expressed uj e M 2 in polar coordinates 

u) = (p cos (p, p simp) e R 2 

and where we note that !F(8 e ) = 9^1r 2 <8><5o- By means of the basic identities cos 2 (<p — 9) + sin 2 (tp — 
9) = 1 and cos(2(y>— 0)) — 2cos 2 (tp—9) — \ we can rewrite operator Bijj in a (second order) Mathieu 
operator (corresponding to the well-known Mathieu equation y"(z) + [(a — 2q) cos(2z)]y(z) = 0, 

MM 

Bu> = D n {{dgf +al- 2qcos{2(<p - 9))) , 

where a = - " +(p2/2 ^ 22+J33) and q = p 2 ( g2 4 2 ^ 33 ) € R. Clearly, this unbounded operator 

(with domain T)(Bu)) = H 2 (T)) is for each fixed u) £ M 2 a symmetric operator of Sturm-Liouville 
type on L 2 (T): 

Bi, = B u . 
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Its right inverse extends to a compact self-adjoint operator on L.2(T) and thereby Bt^ has the 
following complete orthogonal basis of eigen functions 



e%{6)=me n {<p-6,q), neZ,q = p 2 ( 



2 { D 22 ~-Ps3 
4-Dn 



whose eigen-values equal A£ = —a n (q)Du — ^-(I?22 + -D33) < —n 2 Dn < 0, where me„(z,g) = 
ce„(z, q) + ise n (z, q) denotes the well-known Mathieu function (with discrete Floquet exponent 
v = n), |42] . [Tj . and characteristic values a n (q) which are countable solutions of the corresponding 
characteristic equations |12],[T]p.723, containing continued fractions. Note that at u> = 0, i.e. 
p = 0, we have me n (z, 0) = e mz , A° = n 2 . 

The functions g 1— > a n (q) are analytic on the real line. Here we note that in contrast with 
the eigen function decomposition of the generator of the Forward Kolmogorov equation (2.121 of 
Mumford's direction process [19] Green's functions of the contour completion case [19], we have 
gel rather than q € £R and therefor we will not meet any nasty branching points of a n . The 
Taylor expansion of a n (q) for n 7^ 1, 2, 3 (for the cases n = 1,2,3 see pQp.730) is given by 



a„(g) = n 2 + 



5ra 
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9n 2 + 58n 2 + 29 



2(n 2 -l) y 32(n 2 - l) 3 (n 2 -4)* 64(n 2 - l) 5 (n 2 - 4)(n 2 - 9)q e 



0(q 8 ). 



For each fixed u> G R 2 the set {6^f}nez is a complete orthogonal basis for L2(T) and moreover we 
have 



(5 ^) = ^(o)= £ (e^)e^(o) 



for all test functions <j> € P(T). Consequently, the unique solutions of (4.20) are given by 

oo 

K.(u,9)= E W)6^(0)e^, 



n— oo 

71 — — OO 



(4.21) 



Or more explicitly formulated: 

Theorem 4.1. Let Dn, D 22 , -D33 > 0, then the heat kernels 11 D22 ' D33 on the Euclidean motion 
group which satisfy 



{ 9tK = {D n (de) 2 + D 22 {dtf + D 33 (0 V ) 2 ) K 
K.(-,-,0,t) = K,(-,-,2n,t) for all t > 0. 
£(-,-, -,0) =6 e 
/C(-,t) G Li(G), for allt> 0. 



(4.22) 



are given by 



K{x, y, e i9 ,t) = ^ IC^ d ^(lj, e w )](b u b 2 



wht 



Du,D 22 ,D 3 



me n (ip, q)me n (<p - 9, q) taMDll 



2n 



with q = p ^^dJ^ 33 ^ end a n (q) the Mathieu Characteristic (with Floquet exponent n) and with 
the property that /Cf 



-D llt D 22 ,D 33 > Q 



\\K 



T>u,D 22 ,D 3 



lLi(SB(2)) 



= J tf xuD ^ D33 {0,e i6 )&6= C 2 ^)' 1 J ' 



9 d0e~ tn Dl1 = 1. 
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Consider the case where D\\ j 0, then a n (q) ~ — 2q as q — » oo and we have 
lim K. Dll ' D22 ' D33 (u} e lS ) ^ e -i( D 22+D 3 3)(^l+^ 2 y ) e -U D 22-D33)(^l-^ 2 y )§e 

= e -t(D 22 ul+D 33 ^ v ) d e = jc t - D >" D **(u,e ie )6% . 
Finally we notice that the case L*n = yields the following operation on L^G): 

fcO,D 22 ,D 33 ^ se(2) u){g) = J G f^^33 (E -l (x _ x / )) ^ (x / )e ^ ) ^ =( ^ eieG ^ 2 ,^3 , R2 /)(x) 

K 2 

g = (x, e w ) e SE(2) , where Gf 22 ,D33 (x, y) = Gf^ (x) Gf^ (y) equals the well-known anisotropic 
Gaussian kernel or heat-kernel on R™, and where TZ e ie(f>(x.) — 0(i? e " 1 x) is the left regular action of 
50(2) in h 2 (R 2 ), which corresponds to anisotropic diffusion in each fixed orientation layer U(-, -,8) 
where the axes of anisotropy coincide with the £ and fy-axis. This operation is for example used in 
image analysis in the framework of channel smoothing [24], [18]. We stress that also the diffusion 
kernels with D\\ > are interesting for computer vision applications such as the frameworks 
of tensor voting, channel representations and invertible orientation scores as they allow different 
orientation layers {U(-, ■, 0)}eg[o ,27r) to interfere. See Figure [l7j (with inclusion of curvature as we 
will explain in subsection 6.3.11 dependent heat-kernel fcf 11 ' 22,£>33 ; D 22 ^S> D 33 > on SE(2). 
For illustration of the corresponding resolvent kernel R a ,D (with comparison to approximations 
we shall derive in section 4.2 1 see Figure [5] 

Next we shall derive a more suitable expression than (4.21 1 for the resolvent kernel R a ,D- To 
this end we will unwrap the torus to K and replace the periodic boundary condition in 9 by an 
absorbing boundary condition at infinity. Afterwards we shall construct the true periodic solution 
by explicitly computing (using the Floquet theorem) the series consisting of (rapidly decreasing) 
27r-shifts of the solution with absorbing condition at infinity. 

In our explicit formulae for the resolvent kernel R a> D we shall make use of the non-periodic 
complex- valued Mathicu function which is a solution of the Mathieu equation 

y"{z) + [(a - 2q) cos(2z)]y{z) = 0, a,qeR (4.23) 

and which is by definition^} [423p.ll5, [TJp.732, given by 

me±„(z, q) = ce„(z, q) ± ise v (z, q). (4.24) 

Here v = v(a,q) equals the Floquet exponent (due to the Floquet Theorem [42 p. 101) of the 
solution, which means that 

me±„(z + vr, q) = e llJZ mc ±l ,(z, q), (4.25) 

for all z, q £ K. 

Theorem 4.2. Let a > 0, D 22 > D 33 > 0, D n > 0. The solution R^ D : R 3 \ {0, 0, 0} -> K of 
the problem 

(-D n (d e ) 2 - D 22 {d i ) 2 - D 33 {d n f + a) R™ D = aS e , 
R^ D {-, •,#)—*• uniformly on compacta as \9\ — ► oo 



is given by 

R^ D (x ) y,e)=^- 1 [(uj X) oj y )^R^ D (u; x ,uj y ,e)](x,y). 

8 There exist several definitions of Mathieu solutions, for an overview see [l]p.744, Table 20.10 each with different 
normalizations. In this article we always follow the consistent conventions by Meixner and Schaefke | 42| . However, 
for example Mathematica 5.2 chooses an unspecified convention. This requires slight modification of | |4,24[ |, see 
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In case D 33 < D 22 we have 



me, (y, {D ™^ )p2 )mei u (<p - 



(4.26) 



+ me^ („, ^r )p2 ) - ^ (D22 4 ^r )p2 ) uH 



uiii/i u; = (pcos 0, psin 0), where 9 1— > u(#) denotes the unit step function, which is given by u(9) = 
1 i/0 > 0, u(0) = z/6> < and w/iere tte F/ogue£ exponent eguaZs 1/ ( -( a +( 1 / 2 ^ 2 +D 33 ) P 2 ) ^ (d 22 -d^) p 2 ^ 



and whereW a . q — ce„(0, q)se' u (0 : q) equals the Wronskian of ce(- , q) andse(-,q) with a = ( a +( 1 / 2 )CP22+£>33)p 



and q 



_ (D22-D 3 3)P 

4-Di 



In case -D 22 = D 33 (which follows by taking the limit Z? 22 — * D33 in (4--26)) we have 



a e 



which yields for D 33 = D22 '■ 



P = \\u\\,D 2 2 = D 33 , 



R% D (x,9 



-ire 



V On +B 22 



r=||d|,Z)22 = £>33, 



47r v / oTT-D2 



Proof Again we apply Fourier transform with respect to R 2 only, this yields 

0D 22 (%) 2 + Ais^) 2 + i?ii(9 e ) 2 - aI)R™ D = -aS e «■ 

(-I? 22 p 2 cos 2 (<p - 0) - # 3 3P 2 sin 2 ((^ - 0) + Z?n(9 e ) 2 - a/)Ji£ D - -^5° & 

(-D 33 p 2 + (D 33 - £ 22 V 2 cos 2 (^ -9) + D u (d e ) 2 - aI)R™ D = -£6 e o 
((3 e ) 2 + a/ - 2gcos(2(<£ - - -3^*0 



(4.27) 



(4.28) 



where a = - ( a+(p2/2 ^ 2+ ^ ) ) and 9 = p 2 (^f^ 

We shall first deal with the cases -D33 < D 22 and return to the case D 22 = D33 later. In order 
to solve the last equation of ( 4.28[ ), we first find the solutions F, G of the equations 



{{dg) 2 + al - 2gcos(2(0 - 9)))F(9) = 
F(9) -> as 9 -> +00 



and 



((<9 e ) 2 + al~ 2qcos(2(0 - 9)))G{9) = 
G{9) -> as 6> -> -00 



and then we make a continuous (but not differcntiable) fit of these solutions. Now for a < < q 
and a < — 2q we hav^]lm(i/(a, g)) > 0. We indeed have a < and q > since a > 0, £> 22 — D33 > 
and moreover we have 

f< = _ [ ' a+(p 2 /2)(D 22 + D 33 ) \ < _ / (p 2 /2)(Q 22 -D 33 ) \ = _ 2<? 

On / V Dn J 



So consequently (recall (4.25)) we find 



F(9) = Cime-^Op - 0, 9) for 9 > and G(0) = C 2 iae v (<p - 0, g) for 9 < 0, 
now in order to make a continuous fit at = we set C\ — Ame„(^, q) and C 2 = Ame_„(p, g) for 

9 Floquet exponents always exponents come in conjugate pairs, therefor throughout this paper we set the imag- 
inary part of the Floquet-exponent to a positive value. 
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a constant A G R yet to be determined. 



= (9$ - 2gcos(2(^ - 9)) + aI)R™ D Ju>, ■) 

= A(<9 2 - 2qcos(2(ip - 9)) + a/)(me iy (( ( 5, g)me_„(v3 - 0,q)u(6) + me_ u (ip,q)me v (<p - 0,q)u(-0)) 
= -X(me- V (tp,q)mef v (<p,q) - me v (ip, q)me!_ v {<p, q))5 + + + \R™ Dll (5' - S' a ) ^ 

A = -sfe^Kl 1 .?). 1116 -^-.?)])" 1 «■ 

A=-i^(^K(-.9). se -(- ! 9)])" 1 ^ 
A=-ssir(^(0,«)seUO,«)-0) 



(4.29) 



where the Wronskian is given by W[f,g](z) = f(z)g'(z) — g(z)f'(z), which is for solutions of the 
Mathieu-equation independent of z so substitute z = 0. 

Now that we have explicitly derived the solution for the case D22 > D33. We ca n tak e the 
limit D33 t D22 and consequently q J, 0. It directly follows from the Mathieu equation (4.23) that 

limme„(a,?)(0) = me„(a,0) = e 1 ^ 8 for all 9 G [0,2tt). Thereby we have 

9 10 



lim _R a D (u, 9) — — e v 1 1 = 

O33T-D22 ' 47r v A Dn \jD22p 2 + a 



Now the results (4.271 follow by direct computation. 



□ 



We note that if £> 2 2 = A33 the diffusion in the spatial part is isotropic and A = 9| + <9 2 = <9 2 + <9 2 
commutes dg with so in case D22 = -D33 left-invariant diffusion on M 2 x T (with direct product) 
left- invariant diffusion onM 2 xT (with semi-direct product) and the kernels (4--21) indeed coincide 
with the Green's- functions for anisotropic diffusion on M 3 . We have employed this fact in |28j in 
order to generalize fast Gaussian derivatives on images 



/ G L 2 (E 2 ) 



dx m dy n 



(G.*r»/)(x) = (G^* r3 /)(x)gR > 



x G M 2 ,s = (si,s 2 ) G M+xM+, 

(4.30) 



with separable Gaussian kernels Gf =2 (x, y) — Gf^ 1 (x)G~~ 1 (y) (a property which is very useful to 
reduce the computation time) to fast Gaussian derivatives on orientation scores. To this end we 
note that (4.30) can at least formally be written as 



dx m dy n 



e sA f = 



dx r 



d 
dy n 



e sd 'yf 



Now since [A,9g] = we can perform a similar trick for left-invariant Gaussian derivatives on 
orientation scores: 



3 s(Dnd%+D2 2 A)jj 



„sD 22 A 



1_ 

d9 1 ' 



^sDud, 



(4.31) 



d£ m d77"d6» ; " d£ m dn 
which can again be used to reduce computation time: 

_<f^_ {K D 33 =D 22 ^ £(2) tf)( X)e W) = f ( C os9d x + sm9d y ) m (-sm9d x + cosd y ) n G d s ^ 2 22 (x- x')x 



(/-* wG^je - 9>)U{*>.9>)A9>) dx', 



(4.32) 



where we stress that the order of the derivatives matters. 



1G 



Finally, we stress that we can expand the exact Green's function R a ,Dn as an infinite sum 



over 27r-shifts of the solution i?°° n for the unbounded case: 



R a , D {x,y,e ie ) = J™ Y, B£ D (x,y,9-2kn). (4.33) 

N—>oo £ — ' 
k=-N 

Note that this splits the probability-density of finding a random walker (whose traveling time 
is negatively exponentially distributed s ~ NE(a)) in SE(2) at position (regardless its traveling 
time) (x, y, e 18 ) given its starting position and orientation e = (0, 0, e° l ) into the probability density 
of finding a random walker in SE{2) at position (x,y,e l6 ) given it started at e = (0,0, e° l ) and 
given the fact that the homotopy number of its path equals k, for k G Z. 

The nice thing is that the sum in (4.33) (w hich decays rather rapidly) can be computed 



oo 

explicitly by means of the Floquet theorem, i.e. (4.251, and the geometrical series r 

n=0 



for r = e %v with r = \e lv \ < 1 since the imaginary part of v — v(a, q) is positive. By straightforward 
computations this yields the following result. 

Theorem 4.3. Let a,Dn,D 2 2 > and D 33 > 0. Then the solution R a , D : SE(2) — > M of the 
problem 

(-D n (d e ) 2 - A> 2 (^) 2 - D 33 (d v ) 2 + a) R a<D = aS e , 
R a .D{-, + 2kn) = R a , D (; --9) for all k e Z, 
R a . D £ U(SE(2)), 

is given by 

R a , D {x,9) = J2 R aA x > 6 ' + 2fc7 



the righthand side of which can be calculated using Floquet's theorem and (4.26) yielding for 
D 33 < D22: 

[TR a>D {-,0)](u:) = 4^^^(0,9)^(0,,) { 

(— cot(vn) (ce u {ip, q) se„(<p — 6,q) + se„(<p, q) se„(ip — 9, q)) + 

ce„(tp,q) se v {ip - 9,q) - se„(tp,q) ce„(tp - 9,q))u(9) + (4.34) 
(- cot(;/7r) (ce„(</p, q) ce v {ip ~9,q) ~ se„(ip,q) se„(ip - 9, q)) + 

ce„(ip, q) se v {ip - 9,q) + se v ((p, q) ce„(tp - 9, q)) u(-9) } 

with q = ^ D22 4j^ 3 ^ P , <*> = (p cos f, p sin ip) and Floquet exponent v = v{a,q), a = — a+ ( 1 / 2 ) ( ^ 2 - j ~- P33 ) p 
and where 9 1— > u(9) denotes the unit step function, which is given by n(9) = 1 if 9 > 0, u(9) = 
if 9 <0. 

The results in the preceding theory on the resolvent Green's function of the contour enhance- 
ment process can be set in a variational formulation, like the variational formulation in [llj (where 
^33 = 0). 

Corollary 4.4. Let U G L 2 (S'i?(2)) and a, Dn, D22 > 0, D 33 > 0. Then the unique solution of 
the variational problem 



arg mm 

WeB 1 (SE(2)) 

SE(2) 



J ^(W(g)-U(g)) 2 + D 11 (d6W(g)) 2 + D 22 (d i W(g)) 2 + D 3a (dzW(g)) 2 dvsE(2 ) (g) (4.35) 



is given by 

W(g) = {R atD * S E(2)U)(g)= J R a Mh^9)U{h) dp SE(2) {h) 



SE{2) 



where the Green's function R a ,o '■ SE(2) —> R + is explicitly given in Theorem 4-3 
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Proof By convexity of the energy 

S(W) := J ^(W(g)-U(g)) 2 + D n (d e W(g)) 2 + D 22 (d^W(g)) 2 + D 33 (d v W(g))M^ SE{2) (^ 



SE(2) 



the solution of the variational problem (4.351 is unique. Along the minimizer we have 



lim S(W + h5)-S { W) =Q 
HO h 

for all pertubations 5 £ H 1 (SE(2)). So by integration by parts we find 

(a(W - U) - D xl d 2 e W - D 22 d\W - D 33 d 2 W, S) U(SE(2)) = 
for all 5 £ M 1 (SE(2)). Now M 1 (SE(2)) is dense in h 2 (SE(2)) and therefore 

aU=(aI- (D n d 2 g + D 22 d\ + D 33 d 2 )) W 

so W = a (al - (D n d 2 + D 22 d\ + D 33 d 2 )^) U and by left -invariance and linearity this resolvent 

equation is solved by a S'_E(2)-convolution with the smooth Green's function R a: o ■ SE(2)\{e} — * 
R+ from Theorem SH □ 

Remark: We looked for a variational formulation of the contour completion process as well, 
but in vain. 

4.2 The Heisenberg Approximations of the heat-kernels on SE{2) 

If we approximate cos 9 ~ 1 and sin 9 sa 9 the left-invariant vector fields are approximated by 

li = d g , A 2 = dx + 9d y , A 3 = -9d x + d y (4.36) 

which are left-invariant vector fields in a 5 dimensional Nilpotent Lie-algebra of Heisenberg type. 
In our previous related work, we used this replacement to explicitly derive more tangible Green's 
functions which are (surprisingly) good approximation^] of the exact Green's functions of the 
direction process with the goal of contour completion (i.e. a 2 ,Dn ^ 0, other parameters are set 
to zero) for reasonable parameter settings, see [TS]- In fact this replacement will provide Green's 
functions on the group of positions and velocities rather than Green's functions on the group of 
positions and orientations, see [50] App. C. 

Here we will derive the Green's functions for contour enhancement, which are the heat-kernels 
on SE(2). In the case of contour- completion, however, one has the interesting situation that the 
approximative left-invariant vector field A 2 = d x + 9d y together with the diffusion generator {dg) 2 
and the identity operator / and all commutators form an 8-dimensional nil-potent Lie-algebra 
spanned by {I,d x ,dg,dy,9d y ,d 2 ,dgd y ,d 2 }. From this observation and [55] Theorem 3.18.11 p. 243 
it follows that the approximations of the Green's functions (which are again Green's functions but 
of a different Heisenberg type of group of dimension 5) 



pDxi 02=1 - " 3(xe-2y)'+^(G-K i:) 

■ S 



K s lua2 {x,y,9) =6(x~ s) 9n ^ 2 e fa3 ».- 
n°^= ^ _ 3(-.- a „' + -'(.-»o-)' 



where u denotes the lD-Heavy-side/unit step function. This technique can not be applied to the 
diffusion the commutators of the separate diffusion generators provide infinitely many 

directions. Here we follow [llj and apply a coordinates transformation 

= »'•<'> - *• ivm-Tm- tSS) (4 - 38) 



10 In fact in the field of image analysis the approximative Green's functions is often mistaken for the exact Green's 
functions. 
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where we note 



dsK Dll ' D22 — /n a2 1 n ra 1 fla r; " D|l ' D! 



= (Du$ + D 22 (& + 00„r) Kir 11 ' 22 & 



d s K. 



D11.-D22 



Dii,D 22 _ 1 



which provides us the left-invariant evolution equation on the usual Heisenberg group H3 gen- 
erated by Kohn's Laplacian. The Heat-kernel on H(3) is well-known, for explicit derivations see 

and is s iven b y 



2r 



(2tts) 2 J r sinh(2r) \sVD n D 22 



2rt 



g 2tanh(2r) 



dr, 



(4.39) 



and as a result by (4.381 we obtairj^jjhc following Heisenberg-type approximation of the Green's 
function 



K Dll ' D22 (x v 6) = x - K Hs ( x 



8D 
„-l 



9 2(y-^) 

V2D7I' VD 1± D 22 

1 f ?r / *T(y 2 > ) p tanh( 2T ) — 

nD 22 ir 2 s 2 J sinh(2r) COb \s^D 11 D 22 J 6 



dr 



(4.40) 



lim a R a (x, y, I 



i^D 11 D 22 I i , 2 $2 y (y _l^)2 

ie\D 22 + d 11 J "f DnD 22 



See Figure|5]for illustrations of both the exact resolvent Green's function i?^ and its approximation 



4.2.1 The Hormander condition and the Underlying Stochastics of the Heisenberg 
approximation of the Diffusion process on SE(2) 

A differential operator L defined on a manifold M of dimension n € N, n < oo is called hypo- 
elliptic if for all distributions / defined on an open subset of M such that Lf is C°° (smooth), / 
must also be C°° . In his paper |36j . Hormander presented a sufficient and essentially necessary 
condition for an operator of the type 

r 

L = c + X + ^(A^) 2 , r<n 

i=l 

where {Xi} are vector fields on M, to be hypo-elliptic. This condition, which we shall refer to as 
the Hormander condition is that among the set 

{X n ,[X n ,X n }, [X h ,[X j2 ,X j3 )], [X jlt [X h ,[X js , . . .,X jk ]]] ... ./, {0.1... .,r}} (4.41) 

there exist n which are linearly independent at any given point in M. Note that if M is a Lie-group 
and we restrict ourselves to left-invariant vector fields than it is sufficient to check whether the 
vector fields span the tangent space at the unity element. 

If we apply this theorem to the Forward Kolmogorov equation of the direction process than 
we see that the Hormander condition is satisfied since we have M = SE{2) x M+, Xq = —d s — d^, 
Xi = dg and we have 

dim span{-5 s - % dg, [dg, -d s - c%], [dg, [dg, -d s ~ <%]]} = dimspan{<9 s , dg, d 6 , d v } = 4 

and indeed the Green's function of Mumford's direction process is infinitely differentiable on SE(2), 
see |19j . Similarly the Green's function of the resolvent direction process determined by LR = 6 e , 

11 Note that our approximation of the Green's function on the Euclidean motion group does not coincide with 
the formula by Citti in 
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Figure 5: Top Row: Left: isocontour-plots of the marginals of the exact resolvent Green's function of 
the direction process R a 11 ' a2 given in [19], ko = 0, a = A, flu = ^ . 
curves of the marginals of -R^ 11 (given by ( 4.37 1) and i?J 1:L ' a2=1 . 
the approximation R a . The small difference is best seen in the iso-contours close to zero. In the exact 



Right: a comparison of the level 
Dashed lines denote the level sets of 



case oriented photon may loop, whereas in the approximate case oriented particles must move forward. 
Bottom Row, right: A comparison between the exact Green's function of the resolvent diffusion process 
a = jjg, Da = 0.1, D22 = 0.5 in Theorem 4.3 and the approximate Green's function (in dashed lines) of 
the resolvent process with infinite lifetime lim a _,o <x~ 1 (DuAl + D22A2 — al)~ 1 8 e , d 1± = 0.1, d 22 = d.s given 
by (4.401. Bottom left: 3D- view on a stack of iso-contours (top: approximation, bottom exact) viewed 
along ^-direction. 
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with L = — 9j + Du(de) 2 — 7/ is infinitely differentiablc on SE(2)\{e}, for explicit formulae see 
[IT?] . To this end we set M = SE(2) and we note that 



span{d e , [dg,d^], [dg, [dg,d 5 ]}} = span{dg,d ( ,d v } = C(SE{2)). 
However, in the case of the direction process the Heisenberg approximation of the time dependent 



Green's function (4.371 is singular. This is in contrast to its resolvent kernel where the Laplace 



transform takes care of the missing direction direction in the tangent space: 

dim span{-<9 s - (d x + 9d y ),dg, [dg, -d s - (d x + 9d y )], [dg, [dg, -d a - (d x + 0d y )))} 
= dim span{<9 s + d x ,d y , d s + d x + 8d y } = 3 

By the preceding it follows that this deficit does not occur in the Heisenberg approximation 



(4.401 of the pure diffusion (contour completion) case. This can be understood by the Hormander 



condition, since set S = SE(2) x R + then we have 

dim span{9 s , d x + 9d y , dg, [dg, d x + 9d y ] = d y } = 4. 
This puts us to the following question: 

"Can we get physical insight in the induced smoothing in the remaining directions in the diffusion 

processes on SE(2) generated by hypo-elliptic operators which are not elliptic ?" 

Before we provide an affirmative answer to this question we get inspiration from the heat kernel 



on the 3D-Heisenberg group H3, recall (4.391, which is smooth in all directions, despite the fact 
that diffusion is only done in d x + ojdt and d^ — ircVdirection. Here, the induced smoothness in 
t direction, has an elegant stochastic interpretation. As shown in [3T], the underlying stochastic 
process (with the diffusion equation on H3 as the forward Kolmogorov equation) is given by 

Z{s) = X(s) + 1 W(s) = Z + ey/a, e - Af(0, 1) 

T(s) = 2 J WdX - XdW, s > ( 4 ' 42 ) 


so the random variable Z is a Brownian motion in the complex plane and the random variable T(s) 
measures the deviation from a sample path with respect to a straight path Z(s) — Z + s(Z(s) — Z ) 

by means of the stochastic integral T(s) — 2 J WdX — XdWA To this end we note that foi 
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s 1 ► (x(s),lj(s)) S C°°(R + ,]R 2 ) such that the straight-line from X to X(s) followed by the 
inverse path encloses an oriented surface f2 £ M 2 , we have by Stokes' theorem that 2/i(f2) = 
- $*(-X'(t)W{t) + X{t)W'{t)) dt + = Jo WdX - XdW. 

Now by the coordinate transformation in ( |4.38 1 we directly deduce that the underlying stochas- 



tic process of the Heisenberg approximation of the diffusion process on SE(2) is given by 

X(s) + i 0(a) = X(0) + i 6(0) + ^~s{e x + ie e ), where e x ~ Af(0, 2D n ),e e - Af(0, 2D 22 ) 

Y( S ) = x(s) 2 e(s) + 1 / S edx ~ xde = / S e(t) - e(o)dt, 

which provides a better understanding of the "implicit smoothing" (by means of the commutators) 
within the Hormander condition of the Heisenberg approximation of the diffusion process on SE(2). 



5 Modes 

The concept of a completion distribution is well-known in image analysis, see for example |49j . 
|60J, [7], [T5]. The idea is simple: Consider two left-invariant stochastic processes on the Euclidean 
motion group, one with forward convection say its forward Kolmogorov equation is generated 
by A and one with the same stochastic behavior but with backward convection, i.e. its forward 

12 A Brownian motion is a.e. not differentiable in the classical sense, nor does the integral in [ ] 4 . 4 2 [ l make sense 
in classical integration theory. 
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Kolmogorov equation is generated by the adjoint of A* of A. Then we want to compute the 
probability that random walker from both stochastic processes collide. This collision probability 
density is given by 

C (u,v) = (A- aI)~ l U{A* - al^W U,W G L 2 (SE(2)) L 1 (SE(2)) 

where U, W are initial distributions. This collision probability is called a completion field as it 
serves as a model for perceptual organization in the sense that elongated local image fragments 
are completed in a more global coherent structure. These initial distributions can for example 
be obtained from an image by means of a well-posed invertible wavelet transform constructed by 
a reducible representation of the Euclidean motion group as explained in [TB]. Alternatives are 
lifting using the interesting framework of curve indicator random fields [6] or (more ad-hoc) by 
putting a limited set of delta distributions after tresholding some end-point detector or putting 
them simply by hand [BIT] . Here we do not go into detail on how these initial distributions can be 
obtained, but only consider the case U = <5(o,o,0 o ) an< ^ W = ^(xi, 2/1,61)1 x uVi G K- In this case we 



obtain by means of (4.37) the following approximations of the completion fields: 



£a,£Hi,«o,«i( Xi 0) = a 2 ((i-a/)" 1 ^) (x,0) ((A'-a/)- 1 ^,-*) (x,tf) 
= T a:RQ>Dll ^ ,e {x,y,9)T a>KltDll ■- X1 , yi ,e 1 (—x,y, —9), 

with corresponding modes, obtained by solving for 

d y {T a ^ ^ Dll . Q fi {x,y,9)T a ^ 1>Dll ■_ Xlt y lt o 1 (—x,y,—9)} = 
de{T a ^ Ko ^ Dll . Oi g (x,y,9)T aiKl ^ Dll ■- Xl ,y 1 ,e 1 (—x, y, -0)} = 

These modes depend on only on the difference n — K\ but not on Dn nor on a: 



(5.43) 



y{x) = x9 + ^(-22,1 + x 1 (9 - 0i)) + ^(3 Vl + x 1 (9 1 - 20 o )) 
+ (^i)(^- Xl )2 ( a „ x)x 2 

6{x) = 9 + 2^(3 yi + - 20 o )) - 3g(2j/i + aifa - 9 )) 
+ ( -^ 1 x(x - Xl )(-3x 2 + 3x lX - x\) , 



(5.44) 



where x G [0, Xl ] and y(0) = 0, 0(0) = O and 1/(11) = yi, 9( Xl ) = -0! and gl(0) = O and ^(n) 
—^i, see Figure [7] These modes are the unique minimizers of the following variational problem 



(y"{x)- (^-g)^) ^ dx | y (0)=0,v(i 1 )=l«,y'(0)=flb,W / (a!i) = -fli} (5.45) 



argmin{£ (y) 







where c(x) = 2Q{x - - a> and t/'(x) = 0(s) + (Kl ~ Ko a )d(a) with 



d(x) = —2a; (a; — Xi) . The variational problem (5.451, for the case hq = «i is indeed the corre- 



sponding (with arclength replaced by x) approximation of the elastica functional in |43 and indeed 
§f (j/) = for all v G V((0,xi)) if and only if ?/ 4 )(x) = (n\ — kq)c^ 2 \x) under the conditions 
y(O)=O,y(x 1 ) = y 1 ,y'{O) = e o ,y'(x 1 ) = -0 1 . 

We note that because of left-invariance with respect to the 5-dimensional Heisenberg type of 
group we have 

$a,/eo,Bii;x / ,0'( :c J &) = S a , K0 ,D li; e(x -x',y-y' - 9' (x - x'),6 - 9'). 

As a result the approximate completion field (and thereby its mode) is not left-invariant on R 2 xi T 
and thereby its marginal is not Euclidean invariant. As a result the formulas do depend/^] on the 
choice of coordinate system {x, y}. 

However, this problem does not arise for the exact completion field 

cgo , gi ,a,D luKo , Kl = a 2 (( A _ aI yi Sxo ,e ) {(A* - al)- 1 S Xlt6l ) , 



3 If {x,y} is aligned with go the result is different then if it is aligned with gi. 
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since by left-invariance of the generator A we have 

(A — aI)~ 1 S go = (A — aI)~ 1 £ go 5 e = C go (A — al)~ 1 5 e and therefore 

£f/igo,^Si>«.-Dii,reo,Ki — £. £79o,fli,a,-Dn,Ko,Ki £ or & jj /j g ><] T. 

Throughout this paper we shall often use the following convention 

Definition 5.5. A curve s i— > j(s) = (x(s),y(s),e ie ^) in SE(2) is called horizontal iff 6{s) = 
arg(a/(s) + iy'(s)). Then 7 is called the lifted curve in SE(2) of the curve s 1— > x(s) — (x{s),y{s)) 
in R 2 . 

Note that this sets a bijection between horizontal curves in SE{2) and curves in M 2 . 

5.1 Elastica curves, geodesies, modes and zero-crossings of completion 
fields 

In his paper Mumford [13] p. 496 showed that the modes of the direction process are given by 
elastica curves which are by definition curves t 1— > x(f) in R 2 , with length L and prescribed 
boundary conditions 

x(0) = xq,x(Z) = xi with prescribed directions , , 

arg(x'(0) + i y'(0)) = 9 and arg(a:'(L) + i y'{L)) = 9 1 1 > 

which minimize the functional 

£ e (x) = / k 2 (s) + e ds, with e = 4a£)n, with n(s) — ||x(s)||, and s > arclength. 
Jo 

Here he uses the following discrete version (with N steps) of the stochastic process 
' e ie( Sk +A s ) =e i(«WW5 Ek ) ] £fe ~7V(0,a 2 ), D n = \a % 

As = with L ~ NE(a),k — 0,...,N —1, 

where the physical dimension of e k and cr = \/2Dn is [L_EA r GT_H'] _1/ ' 2 , for the definition of the 
mode: 

Definition 5.6. The mode of the direction process with parameters ot,D\\ is a curve in M 2 which 
is the point-wise limit of the maximum likelihood curve of the discrete direction process with given 



boundary conditions (5.4-6) 



In his paper, [43], Mumford states that the probability density of a discrete realization, a 

L_ 

N 



polygon of length L, T = [j x(s;), x(si+i) whose sides have length -4 and the 9i are discrete 



Brownian motion scaled down by As = -fe: 



7 < 4 i = t \j i € {0, . . . ,N - 1}, 



£j independent normal random variables with mean 0, standard deviation a = \/2 Dn, equals 
P(T) = " e *=° 

- E 1 #(^ i )/(2 CT 2 )-aL 

= e i=o v / 
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which converges as N — > oo to 

e -Wlti- 2 (s)+^^ withe = AaDlu £) 11 = i(7 2 ) (5.48) 

from which he deduced that the modes (or maximum likelihood curves) are elastica curves. 

The drawback of this construction is that the definition of the mode is obtained by means of a 
discrete approximation. However, Olaf Wittich brought to our attention that the above construc- 
tion is quite similar to the minimization of the Onsager-Machlup functional J — xgScal^)] 
which under sensible conditions yields the asymptotically most probable path in a Brownian mo- 
tion on a manifold M, |47j . which does not require the jump from the continuous to the discrete 
setting and back. This is a point for future investigation. 

Mumford's observation that the modes of the direction process (with parameters Dn > 
and a > 0) coincide with elastica curves (with e = AaDu) raises the following two challenging 
questions: 

1. Is there a connection between the unique curve determined by the zero crossings of dgC 9o ' 9l ' a,Dl1 
and d v C ga <ei-> a > Dl1 , i.e. the unique curve determined by the intersection of the planes 

{g E SE(2) | d e C 30 ' 91 ' a ' Dll (g) = 0} and {g e SE(2) | d^C 90 ' 91 ^ 11 (g) = 0} 



and the (lifted) modes/elastica curves ? Do they coincide, likewise their respective Heisen- 
berg approximations ( ]5.45 1 and (5.44) 



? 



2. What is the connection between this result and the well-known Onsager-Machlup functional 
which describes the asymptotic probability of a diffusion particle on a complete Ricmannian 
manifold staying in a small ball around a given trajectory, [47 . 

Numerical computations seem to indicate that the unique curve induced by the zero crossings of 
dgC = and d v C = closely approximate the elastica curves ! In fact they even seem to coincide, 
see figure [6] The main problem with mathematically underpinning this numerical observation here 
is that we only have elegant formulae for the exact Green's functions in the Fourier domain. This 
problem did not occur in the case of the Heisenberg approximation. 

In the Heisenberg approximation case the intersection of the planes given by 

dyC 90 ' 91 ' a > Dll {x,y,6) = and d e C 9o > 9l > a ' Dll (x,y,9) = 



yield the i?-spline solution ( |5.44 1 (minimizing the approximate elastica functional 5.45 where the 



role of arc-length s > is replaced by x), which does not depend on a, illustrated in Figure [7] 
This is due to the fact that the approximate resolvent Green's function 

R a ( x ,y,6)=a— ? e ax e ** s °n ux , 

satisfies 

-=•011,02 = 1/ 

lg lua3=1 (x,y,0)=ae- a * lim £*= _^ V '°> (5.49) 

a— »o a 

which coincides with the fact that the random walker in the approximate case is not allowed to 
turn (it should always move forward in ir-direction) . 

Regarding the exact case both the intersecting curve of the planes 

{<? e SE(2) | d e C 90 ' 91 ' a ' Dll (g) = 0} and {g e SE(2) \ d v C 9o ' 9l > a ' Dl1 (g) = 0} (5.50) 

and the elastica curves will depend on a. However, the Green's function^] 11 • a2=1 and thereby 
the intersection of the planes (5.501 , only depend on the quotient (after rescaling position 



14 for explicit formula for the exact resolvent Green functions of the direction proces similar to we refer to 

nsa. 
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variables by x i— > ^==x) whereas the elastica curves only depend on the product e = 4a*Z?n. So 
the curves will certainly not coincide for all parameter settings (a, Dn). Furthermore the resolvent 



Green's function of the exact direction process i?^ 11,a2=1 does not satisfy (5.491. However for the 
special case a — > oo and thereby — > the approximate Green's function converges to the 
exact Green's function, see [T5], [53]. Moreover as a — > oo the arc-length s > of the projection 
of the path of an exact random walker of the direction process on the spatial plane (which is 
differentiable, recall from Figure 3]) tends to the initial direction which is along the x-direction 



and as a result the energy (5.451 (where we recall 9(x) — y'{x)) tends to the elastica energy 
J(k(s)) 2 ds for fixed length curves, since for horizontal curves we have k(s) = #'(s), as a — ► oo. 
So by the different dependence on the parameters Dn and a we can only expect the intersection 
curve of the planes dgC — and d v C = to coincide with the elastica curve in the limiting case 
a — ► oo. 

In image analysis applications we typically have that < « 1 in which case the Heisen- 
berg approximation is a good approximation, as a result for these parameter settings the intersec- 
tions of the planes dgC 9o ' 9l ' a ' Dl1 — and d v C 9o,9l ' a,Dl1 = depend very little on the parameter 
and as a result they turn out to be close approximations of the elastica curves. 

This serves as a theoretical motivation for our curve extraction from completion fields C : 



SE(2) —> M. + of orientation scores 5.43 via the zero crossings of dgC and d v C which is useful for 
detecting noisy elongated structures (such as catheters) in many medical image analysis applica- 
tions. 

Finally, we note that the above observations are quite similar to the result in large deviation 
theory, |59j . where the time integrated unconditional Brownian bridge measure on a manifold 
SE{2) uniformly tends to the geodesies which for a — > oo, see Appendix |B| Here, we stress that the 
Brownian bridge measure on the manifold SE(2) is related to the completion fields of the contour 
enhancement process on SE(2) (with D33 > 0) and in limiting case -D33 — > the corresponding 
geodesies (where we restrict ourselves to horizontal curves) coincide with the horizontal minimizers 
of J a/k 2 + eds. These horizontal curves were also reported by Citti and Sarti as geodesies on 
SE(2) and note that these curves are, in contrast to the elastica curves, coordinate independent 
on the space SE{2) ! 

In respectively Section [5] and Section |7.1| we will derive exact formulae for the curvature of 
geodesies and the curvature of elastica curves which are well-known and we provide a numeric 
comparison between the curves. For reasonable parameter settings these curves turn out to be 
quite close to eachother. The well-known problem with elastica curves is that their curvature 
involves Elliptic functions. Nevertheless, it is still possible to integrate the curvature twice and to 
provide analytic formulae for the curves themselves, |53j . The geodesies however do not suffer from 
this problem, but they cause numerical problems in shooting algorithms, because of singularities 
outside the boundary conditions. 

Therefore we also point to Appendix [A] where we explicitly compute the geodesies. Here we 
will order our results in much more abstract and structured way by means of Pfaffian systems. 
Moreover, by applying the Bryant and Griffiths approach [3] on the Marsden-Weinstein reduction 
for Hamiltonian systems admitting a Lie group of symmetries (developed for elastica curves) to 
the geodesies we are able to get nice analytic formulae, which do not seem to appear in literature, 
for the geodesies which have the advantage that they do not involve special functions. 

Moreover, we refer to Appendix [C] for the computation of snakes/actice shape models in SE(2) 
based on completion fields of orientation scores and elastica curves. See figure [6] 

6 The Underlying Differential Geometry: The Cartan Con- 
nection on the principal fiber bundle Pjj = (SE(2), SE(2)/H, ir, 



The goal is to obtain a connection on SE(2) such that the exponential curves are geodesies and 
compute its curvature and torsion. Moreover we would like to relate the connection to a left- 
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Figure 6: We computed an elastica curve that passes two giv en po ints go = (xo,e 4e °) = (0,0,1) 
and gi = (xi,e 101 ) = (3,-1, j^) via the shooting algorithm (7.95), D\\ = a = j^, and we 
lifted this curve in the Euclidean motion group via 9(s) = Z(x(s), e x ), with parameter e = AaDu 
illustrated by a line of white balls centered around equidistant points on the elastica. Furthermore 
we computed the zero-crossing of the planes where the exact completion field C^ ^^ Ko=Kl=0 = 
a 2 (al + d^ — dg)~ 1 S go (aI ~ c% — dg)~ 1 5 gi (whose xy-marginal we depicted at the bottom) has zero 
respectively 8 (red-plane) and r\ (yellow-plane) derivative. Note that the zero-crossing of these 
two planes is extremely close to the elastica curve. The green arrows reflect ^.(\ gQ and e t\ gi - 
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Figure 7: The shading in these plots denotes the marginal of the analytic completion field 
approximation (5.43) obtained via integration over 9 for x S (0,2) 

torn, and for 6\ = — 15°, 0°, 15° from left to right. The lines drawn on top of these completion 



S (0,2), y G (-0.2,0.8) i.e. 
'o = 0°, 15°, 30° from top to bot- 



fields are the modes ( 5.44| , the optimal connecting lines. 
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invariant Riemannian metric where the left-invariant vector fields {^4i}f =1 = {dg, 9^, d v } serve as 
a moving frame of reference and we want to compute the covariant derivatives of these left-invariant 
vector fields. To this end we recall the general Cartan-connection construction. In contrast to the 
more familiar Levy-Cevita connection this connection does not require a metric. Nevertheless, it 
is possible to relate this connection to a left invariant metric constructed from the Killing form 
on the Lie algebra T e (SE(2)). Unfortunately this killing form is degenerate on SE(2) therefore 
we will embed the Lie-algebra of SE(2) in the Lie-algebra of SO(3) where the Killing form is 
non-degenerate. Throughout this section we will use the Einstein summation convention. 

Let G be a Lie group of finite dimension n with unit element e and subgroup H . By setting the 
equivalence relation a ~ b aT 1 b G H on G we get the left cosets as equivalence classes [g] = gH. 
Let G/H denote the partition of left cosets on G. Let tt : G — > G/H be the projection of G onto 
G/H given by ir(g) = [g] and let R be the right multiplication given by R^g = gh. Note that 
Tr(Rhg) = 7r(ff). This yields a principal fiber bundle Pjj — (G, G/H, tt, R) with structure group H. 
A Cartan/Ehresmann connection^ ui : TP — > T e {G) is a Lie-algebra valued 1-form on Pjj such 
that 

1) Lj((R h )*Y) = Ad(/i- 1 )w(F) for all h G H and all vector fields YonP 

2) w{dK(X)) = X for all X G T e (H), ^ ' 

where (R g )* denotes the push- forward of the right-multiplication and 

Ad(g) = (Rg-iL g ) t (6.52) 

which equals the derivative (at the unity element) of the conjugation automorphism on/n-^ hgh^ 1 , 
h,g G H. Note that requirement 1) means uj g h{(.Rh) *Y g ) = Ad(h~ 1 )oj g (Y g ) , for all vector fields Y 
and all g, h G H. 

In particular we take the Cartan-Maurer form u>h = {L7 )*w e : T g (Pn) — * T e (H) C T e (G), 
,g G G, with uj e = I. By using the restrictions {-4i| g }" = i of the left-invariant vector fields {^}™ =1 , 
with corresponding covectors {d^ l }" =1 with (d^l l ,^j) = Sj (also known as the Maurer Cartan 
co-frame), to g as a local basis for T g (G) for all g G G, where we assume that {^4j}™ =1 are ordered 
such that the first m < n, {Ai, . . . ,A m } elements generate the subgroup H we can express the 
Cartan-Maurer form on Ph as follows 

u g (X g ) = ^(dX J | ff 7 X g ) Tg{G) A u (6.53) 

for all vector fields on Pjj and where we recall that Ai\ g=e = Ai G T e (G). 

Next we give a brief derivation of ( |6.53 1. First recall that the left-invariant vector fields {-4i}™ =1 



satisfy Ai\ g = (L g )*Ai, i.e. they are obtained from T e (G) by push forward of the left multiplication 
and therefor the dual elements (the corresponding co-vector fields) are obtained by the pull-back 
from T e (G) 

&A l \ g = (L g )*dA l 

since we have ((L g ydA l 1 (L g )*Ai) = (dA t 1 Aj) = 5j. Now direct computation yields 

m 

(L h -i).X g =Y l {A i ,{L h -^X g )A i 

i=l 
m 

= E<(^)*dA\A s )^ 

i=l 
m 

i=l 



15 In the common case of Riemannian geometry, with Riemannian connection V, one can create a Lie-algebra 
valued one-form by means of u)(a l d{) = T^dx'(a l di)di c = u)*(a l (9i), where the 1-forms ui^ are given by uj^ = r^-dz 1 , 
where the Christoffel symbols are given by r*j = (dx l , VxjXj), with (dx t ,Xj) = <5* . Necessary and sufficient 
conditions for V the map V given by Vjf^j = ^(X)Xk to be a Riemannian connection are ciP Aw* =0 and 
dgij = gtj^i + 9ik w j- Note however that a Cartan connection in contrast to the Riemannian connection does not 
require a metric. 
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In case of the Maurer-Cartan form we see that 2) is satisfied since left-invariant vector fields 
are obtained by the derivative of the right representation and satisfy X e = (L h -i)*(Lf l )*X e = 
(Lh-i)*Xh = u(Xh). 

Finally we note that in the case of the Maurer-Cartan connection 2) is also satisfied as for all 
h, g £ H and all vector fields Y on Ph we get 

u gh ((R h )*Y g ) = (L {gh) -i)*((R h )*Y g ) = (L h -! o L g -i).((R h ).Y g ) = (L h -i), o (L g -i). o (R h )*Y g 

= (L h -i o R h ) t (L g -i) m Y a = Adih-^ujgYg, 

where we note that R^ and L g commute for any pair of elements g,h £ G. 

The horizontal subspace of T g (G) is defined as TC g = Ker(u; ff ) = span™ =rra+1 {^4 i }. A smooth 
curve within 7 : [0, 1] — > Pjj is horizontal if all tangent vectors c'(t) are horizontal, that is within 
TL c (t)- A horizontal lift c* : [0,1] — > Ph of a curve c : [0,1] — > is a horizontal curve with 
7r(c*) = c. It can be shown that a horizontal lift c* of c is uniquely determined by c* = go and 



7r(po) = Co for some given point go £ G. From the first property (6.511 of the Cartan form it 
follows that Hh g — (Rh)*H g and horizontal lifts arc uniquely determined by right action of H 
in the principal fiber bundle Ph , where T g (G) = V g H g , with V g the space of vertical tangent 
vectors. Consequently, the dimension oiTt g equals the dimension of (G/H). 

Example: G = SE(2) = M 2 x T, H = T and u> = (£(0.0, e^)- 1 )* — (-^(0,0, e-* 9 ))* we have 
V g = span{dfl| g } and TL g = span{9^| ff , drj\ g } and horizontal lifts are obtained by multiplication 
with (0,0,e _ie ) from the right. 

Example: G = SE(2) = R 2 x T, H = Y := {(0,h,e l0 ) ft 6 1} and uj = (L(o,/,,o)-i)*, so 
in components the Cartan-Maurer form reads 

Ug(Xg) = g ,Xg)dy (6.54) 

we have V g = span{9,,| } and Ti g — span{9^| g , de\ g } and horizontal lifts are obtained by multi- 
plication with (0,y, 0) from the right. 

Now that horizontal lifts are determined by the right action of H on G we can introduce the 
concept of parallel transport. To this end we will use the left-invariant vector fields as a frame 
of reference in T(G), i.e. we use their restrictions to g S G, {^4i}" = il g as a basis for T e (G) for 
all g £ G. Now the Parallel transport of a tangent vector X[ ff ] = Ai\ on G/H along a curve 
c : [0, 1] -> G/ff is 

T t (X s ) = c*(£)£\ where c* is a horizontal lift of c and c*(t) = (dA\ ^° ,| ) A- 

This definition is independent on the choice of horizontal lift and r t is an isomorphism between 
the tangent spaces T c iq\{G/H) and T c n\(G/H), Now the covariant derivative of the vector field 
Y on G/H along the curve c : [0, 1] — > G/i? in the point [3] = c(0) is defined as 

Vx,# = lim ^((r'T ^(M - t% 0) ) = lim^rV 1 ^) - Y [g] ), 

with = gf(0). The vector field 1" is called parallel along the curve c if ^ x c{t) Y = 0, with 
X c (t) = c (i) = if * € [0, 1). A curve which is covariantly constant, i.e. 

ViC=0 (6.55) 

is called an auto-parallel. In Ricmannian geometry, such curves coincide with geodesies, i.e. the 
unique smooth curves c, with c(a) = p, c(6) — q which minimize 

f ^(x(c(t)))?(t)iJ(t)dt. (6.56) 
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However, due to the torsion in the Cartan connection auto-parallels and geodesies no longer 
coincide. Moreover, the meaning of geodesies as paths with minimal arc-length can not be used 
as we did not consider a connection induced by a Riemannian metric (yet). 

If we want to express the connection in a (possibly degenerate) metric we note that this metric 



must be left-invariant because of 2) in (6.511. Moreover, because of 1) in (6.51) this metric must 



be invariant under the adjoint representation Ad : G — > T e {G) of G on its Lie-algebra. This brings 
us to the (possibly degenerate) left-invariant metric mx induced by the Killing form K: 



m K (A g ,B g ) = K(A e ,B e ) , for all A, Be C(G) 
K(A, B) = trace (adioadB), 



(6.57) 



where the adjoint representation ad : T e (G) — > B(T e (G)) of the Lie-algebra on itself is given by 
(ad( A))(B) = [—A,B], which is the derivative of the Ad representation mentioned before. For 
the moment we will assume that the killing form is non-degenerate (which is not the case for 
G = SE(2)). The matrix elements with respect to our moving frame of reference AA = &H(Ai)\ 

\g 

the components of mx are given by 



gij = trace((ad„4i o adAj)) 



(dA k ,4 k 4A q ) 



(6.58) 



where the structure constants c\ k are defined by [Ai, Aj] — c k ,A k 



6.1 Vector bundles 

If we consider the trivial case H = e, in which case we have G/H = G, u> — and thereby every 
tangent vector is horizontal, so it does not make a lot of sense to consider a principal fiber bundle 
Ph with structure group H. In such a situation one rather considers the action of the group G 
onto itself. The Cartan form on for example SE(2) would now be given by 



u g (X g ) = (de\ g ,X g )de + (dt\ gl X g )d x 



(dr,\,X g )d v , 



which corresponds to (6.53) the connection in case H = G, but now defined on G rather than 



G/G = {e}. This means that we shall consider the vector bundle E = (G, T(G)), which we shall 
consider next. 

Let t h- > c(t) be a smooth curve in G with X(t) = c n (t)Ai- Let /i : G — > E be a section in E. 
Let {/ifc} be the sections in E aligned with the left-invariant vector fields {.Afc} on G. Then the 
Cartan connection in components reads 



(Dn)X(t) := D x{t) n{c{t)) = a k (c(t))^ k (c(t)) + c l a fc (i)L^( C (t))^-( C (t)), 



(6.59) 



for all sections /i(c(t)) = a k (c(t)) /j, k (c(t)) and where D^./ij = Y k ^ k . The Cartan connection D is 
given by D = d + ui since 

DaVfc = d(a k )fi k + a k u(fi k ), 

where we note that by the chain rule we have 

da k 



(d(a k )^ k )X(t) = ^-^ k (c(t))c l (t)S k = ^ k (c(t))^c k (t) = ^ k (c(t))a k (t), 

where we used short notation a(t) = a(c(t)) and moreover 

(a^ic'Ai) = a k (t)^(c(t)M^) = t i (t)a k (t)n j (c(t))Tj k (c(t)) 



(6.60) 



(6.61) 



and the terms in (6.60) and (6.61) indeed add up to the right hand side of (6.591. Now the Cartan 



connection can be split up in a symmetric and anti-symmetric part: 

Ki 



' /,/ — Ik) > Kill ~ 2^kl ^lk) 
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where L fc; are the components of a Levy-Cevita connection induced by the metric mx given by 
(6.571. Now by left invariance we have with respect to our frame of reference that 

Tfcz = \g lm {9mk,i + 9mi,k — 9ki,m) = , where we use the convention: indices after the comma to 
index left-invariant differentiation (i.e. g m k,i — Aig m k)- So in our case the components of the 
Cartan tensor coincide with the components of the contorsion tensor K: 

Tfc; = K l kl = l -g m {c mk i + c m .ik - c Mm ) (6.62) 

where c mk i — gmp^f-n where (? kl are the structure constants of the Lie-algebra 

= (? kl A p ,p,k,l = 1, ...,n, 
where the components of the left-invariant metric tensor are given by 

1 , , 1 



g iS (g) =m K (Ai\ g , A 3 \ ) = -^c* jk c& = ^ ik c\. (6-63) 

1 ad * 



where the Dynkin index J a( j of the adjoint representation, coincides with the dual Coxeter number 
[27], which in case of SO (3) = SU{2) = A 1 equals 2, [35]. 
Now the curvature of the Cartan connection D is 

D 2 = D o D = (d + cu) o (d + lo) = duj + lo A u, (6.64) 

where we note that since u> is a one-form we have Z? 2 /i = dw/i — wd/i + wd/i + lo A lo[i for all section 
/i in E. In components (for details see [37Jp. 111-112) (6.641 reads 

D 2 = hwj ,i - Wi, j + [wi,Wj])dC A d£_ j 
This provides the following formula for the curvature tensor: 

nj _ yj _ yj , V" 1 " yj yX _ yj y\ 

— 1 li, k 1 hi, I ' l^\=\ 1 fcA 1 li 1 (A hi 

E n yj y\ _ yj y\ 
\=l L k\ L li MA A fc»> 

which after some computation using the symmetries of the curvature tensor gives 

4 M = (6-65) 

which is the formula by [3] p. 187. 

We would like to apply ( |6.65 1 and (6.58) to the case of the Euclidean motion group. But here a 
problem arises as the metric ttik induced by the killingform K on SE(2) is degenerate. Therefore 
we embed, by means of a small parameter /3, the Lie-algebra of SE(2) spanned by {ds,d^,d v } 
into the Lie-algebra of 50(3) which is so(3) = {X £ M 3x3 | X T = -X} whose Killing form is 
non-degenerate. We simply obtain the appropriate sectional curvatures and covariant derivatives 
by taking the limit 8 — > afterwards. 

The Euler angle parametrization of 50(3) is given by R ez ~jR e gRe z ,a- Then a basis of left- 
invariant vector fields on 50(3) (in Euler-angles) is given by 

B\ = cot 8 cos 7 & — ^2 da + sin 7 ds 

' ' sin p P 

B 2 = - cot ~B sin 7 % - jg[ d & + sin 7 (6.66) 
63 = ftp 

where we note that [Si,Ba] = S3, [82,83] = B\, [S 3 ,Si] = B 2 - Now apply the coordinate trans- 
formation 

ol — 8x 

~B=\- arctan(/%) (6.67) 
7= 6 
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Figure 8: The embedding of SE{2) in SO (3) by ( |6J7| for j3 = 1. The group 50(3) can be 
identified with a unit ball Bo ^tt with radius 27r by means of the Euler angle parametrization. Here 
all points on the sphere {x € K 3 | ||x|| = a = 2tt} are identified with the origin. Now given a fixed 
member of (x, y, e 10 ) € SE(2) we can obtain the corresponding element in 5*0(3) as follows. First 
consider the point (x,y,z) = (x,y,§) with attached direction 9 — Z((x, {/), (1, 0)) and construct 
the unique half- line I trough the origin with direction 9. Then project (%, y, 0) on the y-axis and 
rotate the point so that it ends up at point P at the line I. Then find the unique point on the 
unit sphere Q such that OQJ-PN where N is the north-pole. Finally, scale Q with x modulo 27r. 



and multiply B\ and Bi with /3 then we obtain the following vector fields 



A 2 



= = -(3 2 y cos 6 dg + cos 6 yfi + f3 2 y 2 d x + sin 9 (1 + [3 2 y 2 ) d y 
= f3B 2 =P 2 y sin d e - sin 9 y/l + (5 2 y 2 d x + cos 9 (1 + f3 2 y 2 ) d y 
= B 3 - 8g. 



(6.68) 



These vector fields again form a three dimensional Lie algebra: 

\A P A 13 ] - A 13 \A S) A (i ] - -A p \A P A p ] - B 2 A () 

and which converges to {A\, A%, A3} — {dg, d%, d n } for (3 — > 0. To get a geometrical understanding 
of the embedding of SE{2) into SE(3) by means of the coordinate transformation (6.671 see Figure 
[8] The components of the left-invariant metric tensor are given by 



9ij(g) = m K (Ai\ , Aj\ ) = -trace(ad(Ai) o ad(Aj)) 



(6.69) 



where c l j k denote the structure constants of the perturbed Lie algebra, which are related to the 
structure constants of the Lie algebra so(3) of 50(3) properly scaled by f3. As a result the 
components of the left-invariant Killing-form-metric with respect to this perturbed Lie-algebra, 
recall (6.631, are given by^ 



G = Iga] = 




(6.70) 



3 Note that the physical dimension of e is [length]' 
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and the non-zero components of the constant Riemann-curvature tensor, recall (6.651, are now 
given by 

^212 — ~^221 — — R331 = R313 = R332 = ~^?323 ~ "-^232 = -^223 = \ >0aS/3^0 

E>2 _ p2 _ p3 _ r?3 _ 1 
^112 — — ^121 — ^113 — ~~ -"431 — I 

From which we deduce that in the limiting case (3 — > the non-normalized sectional curvature 
\dg A df:\K(de A d 6 ) = (R(d0,d^)d^,d e ) = R\i\i = gi m R™i 2 of the plane spanned by {d^,dg} 
is constant = \. Similarly the non-normalized curvature of the plane spanned by {djjjdg} is 
constant = ~ in contrast to the spatial plane spanned by {0^,8^}, which is of course flat. This 
explains the presence of curvature in orientation scores. Moreover, the curvedness of the Cartan 
connection on the space SE(2) is important for application of differential geometrical operators 
on orientation scores. The covariant derivative of a vector field v on SE(2) is a (l,l)-tensor field 
whose components are given by 

The covariant derivative of a co-vector field a on SE(2) is (0,2)-tensor field with components: 

Vjdi = djdi - r^ofc. 

3 

In particular for the gradient dU — U^dA 1 of an orientation score U : SE(2) — ► C and the 

»=i 

corresponding vector field 

dU „ dU „ dU 



This yields the following covariant second order derivation:^ 7 ] on orientation scores U (in the 
limiting case e — > 0): 

/ dedeU d s d e U+ldr,U dr,d e U-±d ( U\ 
[ViVjIT] = [V^] = [9,(7,, - r& U, x ] = ded ( U - 17« 

/ I/flfl §(t/fl« + (7« 9 ) \{Ue n + U n e) \ 

= Hu e9 + u,t) u (( u (v 

\ l(U ne + Ue v ) U ni U m J 

(6.71) 



with Uij = dj(diU), i,j = 1, ... ,3, where we recall (6.621 and d v = [de,d^] and 9j = —[dg,d v ]. 

So by anti-symmetry of the Christoffel symbols r* fc = — Tjj. • we have 

Vi S U + VjiU = AiAjU + A 3 AU. (6.72) 

and therefore our linear and non-linear (that is the conductivity (Dij(W))(g, s) explicitly depends 
on W) evolution equations on SE(2) can be straightforwardly expressed in covariant derivatives: 

d s W(g,s) = g Ai(Dij(W))(g, s)AjW = g V^WX*?, ajVj-W 

where |t//(<7, s)| denotes the absolute value of the orientation score Uf G L2(5£'(2)) of im- 
age / G L 2 (K 2 ) where we note that V, . V,i /),,V,r i = V, , D.XXJ • V,,i V./^,, iV,/ = 



17 Note that a rescaling of 9, say 6 — > A0 directly results in a rescaling of the covariant derivatives: VgVj 
AV e V ? , V S V, ; -» AV e V^, V s V e -» A 2 V e V e . 
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6.2 Auto-parallels 

Recall from Ricmannian differential geometry that curves which are covariantly constant (or auto 
parallel) that is 

V k x = <S> x l = -V kl x k x\ 

where x = x l di and where the well-known Christoffel symbols read rjy = \g lm {gmk.i+9mi.k—gki,m) 
coincide with the path-length minimizers, i.e. geodesies, if and only if the connection is torsion 
free. 

The Cartan connection, however, is not torsion free. Therefore the auto-parallels on SE(2) do 
not coincide with the geodesies. In fact the auto-parallels coincide with the exponential curves. 
To this end we note that the Christoffel symbols (6.621 with respect to our basis of left-invariant 
vector fields are anti-symmetric and as a result we have for auto parallel curves 

d 

ds 

and thereby we have that 

(dA k , c'(s)) = (dA k ,c / (0)) = constant = c k for k = 1, . . . , n, 

so c'(s) = c k -4ft| c ( s ), for some constants Ck € M. Now these curves exactly coincide with the 
exponential curves c(s) = exp(sc fc -4fc| c (o)) c (0) within the Lie group. 

Example : 



Recall that the left- invariant vector fields in SE(2) were given by (2.5 1 as a result auto-parallels 
7 in the Euclidean motion group are given by the following set of equations 

A 7 ( s) = c fc A| 7(s) 7(s),fc= 1,2,3, 

or explicitly in (x,y, 9) -coordinates (not to be mistaken with the (£, 77, #)-coordinates) 
73 (s) = 1 

71 (s) = c 2 cos7i(s) - c 3 sin7i(s) 7(0) = g 
72 (s) = c 2 cos 71 (s) + c 3 sin7i(s) 

the unique solution of which is given by 

7(t) = exp(i(X) cMi))5o = (»o + %r(cos(cH + 9 ) - cos6» ) + jr(sin(cit+6lo) - sin6» ), , . 

i=l (D.Mj 

t/o + #(sin(c 1 i+^)-sin0 o )-^(cos(c 1 t+0 o )-cos0o),e^ cl ^ )), 

// c 2\2_|_( c 3\2 

for c 7^ 0, which is a circular spiral with radius — — and central point 

c c c 
( — t cos 6>o 7 sin 9 + xq, —r cos 6 r sin O + ?/o) ■ 

This result is easily deduced by the method of characteristics for first order PDE's. For c 1 — we 
get a straight line in the plane 9 = 9 : 

j(t) = (x + * c 2 cos 6»o - tc 3 sin 6» , yo + tc 2 sm9 + t c 3 cos6» , e* e °), 



which coincides with (6.731 by taking the limit c 1 — > 0. 
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6.3 Fiber bundles and the concept of horizontal curves 

Recall Definition |5.5| where we provide the definition of a horizontal curve in SE(2). Next we will 
set this definition in a differential geometrical context (which justifies the word horizontal). 

In case G = SE{2) = R 2 X T, Y = {(0,h,e M ) del} and to = (L (0)I/)0 )-i)» we have 
V g = span{d r) \ g } and TL g = span{9j| g , dg\ g } and horizontal lifts are obtained by multiplication 
with (0, h, 0) from the right , where we note that 

5(0, h,Q) = (x,y,e ie )(0,h,0) = (x - h sin 9, y + h cos9,e ld ) = g + he, r 

This particular choice is important in image analysis as it is the differential geometrical description 
of "lifting" of curves. That is to each smooth planar curve in C(E + ,K 2 ) we can create a curve in 
C(R+,SE(2)) by setting 



s 1 ^ x(s) G C(M+,R 2 ) 



<— > 



s 1 — ► (x(s), e i0 &) G C(SE(2), E) with 0(s) = arg(x'(s) + iy'(s)) = Zx'(s). ^ 6 - 74 ^ 
Convention: By Za/(s) we mean the angle that a/(s) makes with the fixed e x -axis, i.e. 

a/ (s) • e-s 



0(s) = Za/(s) ^> cos 0(s) 



In this setting such a curve is horizontal, since its tangent field is a horizontal vector field. Notice 
that right multiplication with a fixed element (0, h, 0) provides a horizontal lift of the curve: 

g(s) is horizontal =>- 5(0, h, 0)(s) = <7(s)(0, h, 0) is horizontal . 
7r( fl (a))=7r(ff(s)(0,fc,0)), for all s > 0. 

where we note Zb'(t) = 0(i) => Z(/i0'(i)(- cos - sin 0(i)) + b'(t)) = 0(*). We note that 
in general right multiplication of a horizontal curve with a constant element in SE{2) does not 
preserve the horizontality. Actually this only works for the subgroup Y. This is in contrast with 
left multiplication with a fixed clement: 

g(s) is horizontal =>■ (gg)(s) — g g(s) is horizontal for all g € SE(2). 

The setting in this example is important to relate elastica curves to geodesies which minimize 

d SE{ 2)(g, 3o) = inf vWF+e|RWF dt \ 

7 is a smooth horizontal curve connecting g and go , 7(0) = (70,7(1) = g} , 

(6.75) 

since only for horizontal curves we have «(s) = 0'(s), so we may write 

dsE(2)(g, go) = inf{ / VT^MF+~e ds I 7(0) = 5o, l(L) = g, 7 smooth and horizontal }. 
Jo 

The auto-parallels in the fiber bundle are now given by 

2 c 2 c 2 1 

7(t) =exp(i(^c l A)) 5 o = (a;o + T (sm(c 1 t+6'o)- sin o ),y OT (cos(c 1 t+0 o )-cos0o), e l < c 



i=l 



(6.76) 



where we note that they follow from the vector bundle case by omitting the vertical direction d v , 
so we get them from (6.73 1 by setting C3 = 0. Note that these auto-parallels are indeed horizontal 
as we have 

arg{7i(i) + i7 2 (t)} = arg{c 2 cos(cii + 8 ) + i sin(cii + O )} = c\t + 9 . 
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Figure 9: All horizontal exponential curves through a fixed point g € SE(2) for different cur- 
vature values, shown from 2 different perspectives. The left-sided image shows that these curves 
correspond to circular arcs if projected onto te spatial plane. 



Also the covariant derivatives are again blind for the vertical direction and they are given by 

X7TT] _f dgdeU d £ d e U \_f d e d e U dgd^U \ 
LVjViUj - ^ j-y g ^ sU j 

for horizontal gradients dU = dgl/dO + d^Udt; , i.e. d v U = 0. 



6.3.1 Horizontality and the extraction of spatial curvature from orientation scores 



Orientation scores U and their absolute value \U\ = \Uf\ = yJ($t(Uf)) 2 + (3(J7/)) 2 in general do 
not satisfy d v U = 0, 9^|J7|. Nevertheless, in our linear and non-linear diffusion schemes (in section 
[3] and in section [8]) on orientation scores, we include the direction + ndg, where k equals the 
horizontal curvature (i.e. the spatial curvature of projected curves on the spatial plane). 

Let U : SE{2) — > R + be some positive smooth function on SE(2). This could for example 
be the absolute value of a (processed) orientation score of an image, which is positive and phase 
invariant see Figure [l] (d) . 

Then by embedding SE(2) into R 3 , the exponential curves through go with direction d AiU\ go 
form "tangent spirals" to the orientation score U : SE(2) — > M + . In particular, horizontal expo- 
nential curves are exp(s(ndg + d^))go and given by (6.761. See Figure [9] In this section we will 
obtain fast algorithms for curvature estimation at position, say go G SE(2), in the domain of U, 
by finding the tangent spiral (exponential curve) through go that fits U in an optimal way. For 
the exact definition of such an optimally fitting (horizontal) tangent spiral we first need a few 
preliminaries. 

We introduce the following (left-invariant) first fundamental form on T(SE(2)) x T(SE(2) 



gij dA l <g> dA j = d6®d6 + /3 2 d£ <g> d£ + /3 2 dr] <g> dr] 



(6.77) 



where f3 > 0, gij is the diagonal matrix with {1, f3 2 , f3 2 } as respective diagonal elements. To this end 
we recall recall (6.70), where we note that (6.771 and (6.701 coincide. Recall that the metric (6.771 



does not coincide with the Cartan connection on SE(2) since it is not right-invariant. However, 
the corresponding metric connection does correspond to the Cartan connection on 5*0(3), recall 



(6.661 and (6.681 



The physical dimension of (3 equals 1/length and j3 is the fundamental parameter that relates 
distance on the torus to the distance in the spatial plane. The inner-product between two left- 
invariant vector fields is now given by 
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where we use the convention c\ = c\.,c\ = c|, c| = c2, fc = 1,2. The norm of a left-invariant vector 
field c l „4i is now given by 



|cVl^ = yJ&Ai^Aih = y/(i*)2 + (l3c£) 2 + (l3ci)* =: Hcll/j, 



(6.78) 



with c = (c , c , <r) G 



Here we stress that the norm 



C(SE(2)) R+ is defined on the 



space L(SE(2)) of left-invariant vector fields on SE{2), whereas the norm || • ||^ 
defined on R 3 . 

The gradient dU of f7 : S£(2) ->■ M+ is given by 

Jrr 9Z7,„ <9C/ at/ 

df/ = w d0+ w de+ ^ d?? - 

It is a co-vector field. The corresponding vector field equals 



(6.79) 



where : T(SE(2))* — > T(SE(2)) the inverse of the fundamental bijection between the tangent 
space and its dual. Note that G~ 1 dAk = Q Ai, with g 1 ^ gti = &\&v The norm of a co-vector field 
is given by 

|aiCL4*|| = g ij <H<ij = (a e ) 2 + /3~ 2 (a c ) 2 + /^(a,,) 2 = ||a||^-i with a = (a 1 , a 2 , a 3 ). 

Finally we stress that if we differentiate a smooth function U : SE(2) — > K + along an exponential 
curve 7(f) = <?oe x p(^(E c 1 ^)) passing go we get (by application of the chain rule) 



|f/( 7 (i)) = (dC/, f (*)) = Eli c l A t U\ l(t) 

= c 1 U e ( 1 (t)) + <?U e (j(t)) + c )i U r ,( 1 (t)). 



(6.80) 



Or in words: The exponential curves {g®e tc Ai } go eSE(2) « r e the characteristics of the left-invariant 
vector field c l Ai- 

After these two preliminaries we return to our goal of finding the optimal tangent spiral at 
position .go G SE(2) given U : SE(2) -> R+ . 

Definition 6.7. The solution of the following minimization problem 



c» = are mm 



dU(y(t)) 



r (t) = goexp(t(J2 cM<)) ; (c 6 ) 2 + /3 2 (c f ) 2 + /?V) 2 = 1 , (6.81) 



yields the optimal tangent spiral {go£ tc * } go eSE(2) °^ position g$ G SE{2) given U : SE{2) 



By means of (6.80 1 and the chain rule the energy in (6.81 ) can be rewritten as 



|(dC0(7(*))U o |! =||V(VC/) r ( 7 (0))-7 / (0)||L 1 





d 6 (dgU) 


d v (d e U) \ 




2 




d £ (d 6 U) 


d n (d S U) 


(1) 




d e {d n U) 


dz(d n U) 


d v (dr,U) J 


go \ ' 


/3-i 



l 9 o Cllj- 



(6.82) 

where VJ7 := (dgU, d^U, d^U) and where the non-covariant Hessian HU is not to be mistaken with 
the covariant Hessian form consisting of covariant derivatives of the Cartan connection (6.71 1. We 
return to this later. 



Note that the minimization problem (6.81 1 can now be rewritten as 
argmin{||(i/[/)(. go )c|| 2 _ 1 | ||c||^ = l} . 
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Set Mp := diag^,/? -1 ,/? -1 } e GL(3,R) and ifgi/ = MpHUMp, then by the Euler-Lagrange 



theory the gradient of V c || (JJJ7)c||g_ 1 



V c (c,(HU) T M%(HU)c) 1 at the optimum is linearly 



dependent on the gradient of the side condition V c (l — ||c| 



V c (c,M- 2 c) i: 



(HU(g )) T M 2 (HU(g ))c* = A M^ 2 c* & (H pUf (H pU)c = Ac, 
for some Lagrange multiplier A 6 1, where c = Mr'c,. 



So we have shown that the minimization problem (6.81 1 requires eigensystem anal ysis of 
(HpU) T HpU rather than eigensystem analysis of the covariant Hessian given by (6.711. The 



eigensystem of the covariant Hessian, however, correspond to the Euler-Lagrange equation for the 
following minimization problem (for simplicity we set (3=1) 



arg mm 



dt 2 



£7( 7 (t)) 



7 (t) - <7 exp(t£>%)) ; J>*) 2 = (c 9 ) 2 + (c«) 2 + (c") 2 = 1 j , 

(6.83) 



i=l 



which by means of (6.801 and again the chain rule can be rewritten as 



|rtf(7(*)) =\&(VU. V(f))f = c T (HU)c 



and as a result the Euler-Lagrange equations for the minimization problem (6.831 correspond to 
the eigensystem of \(HU + (HU) T ), which coincides with covariant Hessian VV T C/ given by 



(6.711: 



VV T [/c = -(HU + (HU) T )c = Ac. 
Experiments on images consisting of lines with ground truth curvatures show that minimization 



problem (6.831 is certainly not preferable over (6.811 for spatial curvature estimation 



Remarks : 

• On the commutative group M 2 (i.e. the domain of images / rather than the domain of the ori- 

fxx fxy 



entation scores Uf) we do not have this difference, since here the Hessian Hf ■ 



fyx fy 



is square symmetric and thereby Hf = \(H f+(H f) T ) and (H f) T (H f) have the same eigen- 
vectors with respective eigenvalues {A„} and {(A„) 2 }. 



If the spatial gradient vanishes at go then 9^U\ = d v U\ go = 0, the problems (6.81) and 
(6.831 have the same minimizer and in this case the covariant Hessian (6.71 1 and the non- 
covariant Hessian (6.82) coincide. 



Sofar we did not include the concept of horizontality. Formally, because of the shape of our ad- 
missible vectors/distributions ip in the wavelet transforms, the orientation scores Uf = Wipf and 
their absolute value | Ut \ usually do not have a horizontal gradient at locations go of elongated 
structures, i.e. in general the gradient does not satisfy c^lt//!^ = 0. Nevertheless, our algo- 
rithm in section [8] requires horizontal curvature estimates from the absolute value of a (processed) 
orientation score \U\. 

Therefor we suggest the following 2 methods for curvature estimation: 

1. Compute the eigen vectors of (Hp or \U\) T (Ha or \U\) with horizontal Hessian 

/ (3 2 d e d e \U\ 0d e d e \U\ \ 
H$T\U\ - (3d e d $ \U\ d&\U\ (6.84) 
V pd e d v \U\ d^d v \U\ ) 
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Figure 10: Curvature estimation test images. 



to this end we note/recall that the optimum c* = arg min{||i^ or |L T |(<7o) c ll|-i I || c ll/3 = 1} with 
c = (cV«) = c e e e + c«e<; satisfies 2(H^° r \U\) T H p lor \U\ c = 2Ac, c* = Mpc for some Lagrange 
multiplier A. Then we compute the curvature of the projection x(s(t)) = Pr2 (g exp(t(^ c\ Ai))) 
of the exponential curve in SE(2) on the ground plane from the eigenvector c* = (c*,cf) with 
smallest eigen value: 

n est = ||x(s)||sign(x(s) • e J( ) = -| (6.85) 

c; 

2. An alternative approach, however, would be to compute the best exponential curve where 
we do not restrict ourselves to horizontal curves. In this case we compute the curvature of the 
projection x(s(t)) = Pr2 (goexp{t(Y^H=i c\ Ai))j of the optimal exponential curve in SE(2) on the 

ground plane from an eigenvector c* = (c*,c*,c2). This eigen vector of (Hp\U\) T (Hp\U\), where 
the 3 x 3-Hessian is given by 

/ (3 2 d g dg\U\ 0dtd g \U\ (3d v d e \U\ \ 
H P \U\ = Pdedtpl d&\U\ d v 8 ( \U\ , (6.86) 
V Pdodr,\U\ d t d^\U\ d v d v \U\ J 

belongs to the pair of eigen vectors closest to the plane { e {| go > e e| go } an d has the smallest eigen 
value. The curvature estimation is now given by 

Kest = ||x( S )||sign(x( S ) • e ?) ) - j C * Slgn(C * ) ■ (6.87) 



V(cS) 2 + (c2) 2 



Note that in this alternative approach, in contrast to the other, we do not include the concept of 
horizontality by restricting ourselves to fitting only horizontal exponential curves, but we simply 
discard the eigen value (which may be small) corresponding to the eigen vector which is most 
pointing out the "correct" horizontal plane in our selection of eigen vector with smallest eigen 
value. 

3. In stead of the Hessian in approach 2. one can also use the eigen vectors of the so-called "struc- 
ture tensor", given by G t * (d g \U\, Pd 6 \U\, pdr,\U\) T (d e \U\, f3d € \U\, Pd v \U\), this corresponds to 
the method proposed by van Ginkcl [32 who considered curvature estimation from non-invertible 
orientation scores. 

For curvature estimation (comparing all three above methods) on orientation scores of noisy 



example images, see Figure 10 and Figure 11 



7 Elastica 

Let e > 0. Then for a smooth curve 1 1— > x(i) in K 2 , with length L we define 

S £ (x) = [ n 2 (s) +eds. 
Jo 
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Image (a) 



l/x est Full Hessian l/K es Hessian Horizontal 1 /^ est Structure Tensor 




Image (a)+noise 



l/x est Full Hessian l^estiessian Horizontal 1 /-*est Structure Tensor 




Image (b) 



l/x est Full Hessian 1 /^esHessian Horizontal i/JQst Structure Tensor 




Error Full Hessi 



10 20 30 40 



Horizontal 



Errorstructure Tensor 
0.7 



10 20 30 40 



Figure 11: Top: Curvature estimation results (left column approach 1. 
column approach 2. in section 



in section 



6.3.1 middle 



6.3.1 right column approach 3. (by van Ginkel 32J) in section 



6.3.11 on image (a) and image (a)+noise for the three methods. For both images, the first row 



shows the density plot of true curvature against estimated curvature, and the second row shows the 
^2-error as function of the different curvatures. Bottom row: Shows curvature estimation results 
on image (b) (see Figure 10 1 
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where s denotes the arc-length parameter^ defined by 



s(t) 



Idr, 



and where curvature of the planar curve is given by k(s) = ||x(s)||. 

Let C = {s ^ x(s) € C°°(IR + ,IR 2 ) | x(0) = x ,x(l) = x 1; (x = 6 , (x = X } be the space of 
smooth planar curves which connect x and Xi such that the starting and ending direction are 
prescribed. 

We sometimes also consider Cl = {s i— > x(s) G C | s G [0, L]} the space of smooth planar curves 
with /ta;e<i total length L which connect Xo and xi such that the starting and ending direction are 
prescribed. 

On C we can consider the following optimization problem, for e > 0, 



Find s i — ► x(s) £ C such that £ e (x) is minimal. 
Similarly we can consider the optimization problem on Cl- 
Find s i— > x(s) S Cl such that £o( x ) is minimal. 



(7. 



(7.89) 



We first consider (7:. 



so let e > 0. Let x be the optimal curve with tangent t = x, normal 
n = x and curvature n = x ■ n. Then any infinitesimal deformation of this curve should yield lower 
energy. Since we can always re-parameterize our curves we only need to consider deformation of 
the curve in normal direction 

xnew(s) = x(s) + hS(s) n(s),h > (7.90) 

with 5 twice differentiable and compactly supported within the open interval (0,L). Then we 
stress that the arc-length parameter s new of the pertubed curve xjv ew does not coincide with 
the arc-length parameter s of the original curve x. In fact we have 

x New (s) = (1 - ShK(s))(t(s) + hS'(s) n(s)) + 0(h 2 ) = (1 - 5hn(s))(t(s) + hS'(s) n(s)) 



Then 



ds 7v ew = (1 — hS n)ds 



{•NEW 

knew 



dSNEW 
d Xjy euj 



ds n ew ds 

n - hd't 



hS'n 



dt 



knew = 



NEW 



ds 



NEW 



to-new = n+hS"+hSK 2 +0(h 2 ) 



(7.91) 



From which it follows that 

hio h 



hfo ^ (■l'o NEW k new( s ) + dsNEW - Jo k 2 (s) + eds) 
Km I J n L ((n + hS" + hSn 2 ) 2 + e) (1 - hSn) - {n 2 + e)ds 



so by partial integration we find lim 



£{x+hS n)-£(x+5 n) _ 
h 



= for all 5 if and only if 



2k- 



(7.92) 
They can 



however we stress that not all solutions of (7.921 lead to global minimization of (7. 
be local minima or even saddle points. 

For problem (7.891 we note that the deformations must be length preserving, in this case we 
have 



<1-Snew = / ds <t=> / n(s)5(s)ds = 



3 arc-length in M 2 , not arc-length in SE(2). 
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Consequently the optimation is similar as above with the only difference that e has to be replaced 
by an Euler lagrange multiplier 



2k 



= Xk, 



A > 0. 



(7.93) 



Now (7.921 and ( |7.93[ ) provide the curvature of elastica curves, which is unique if we set 

k(0) = kq and k(0) = K , 

for some positive constants ko and k' which we will determine later. To get the elastica curves 
themselves we have to integrate the Frenet formulas: 



d_ 
ds 



x( S ) 
x( S ) 





-k(s)J 2 



k(s)I 2 




x( s ) 
x( S ) 



with I2 the identity matrix and where the solution is uniquely determined by s 1— > k(s) and 
Zx(0) = 9q, to this end we note that the exponential of a skew symmetric matrix is orthogonal 
and therefor ||x(s)|| = ||x(s)|| = 1 and x ■ x = for all s > 0, so that Zx(0) = 9q sets the initial 
condition x(0) = cos^oe^ + sin 0Qe y , x(0) = — sin#o e iE + cos#o e y and thereby the full solution 
s i-> (x(s),x(s)). 

Now we get the elastica curve by integration over s > 0: 

x(s) = x(0) + J x(t) dt. 


The three free parameters k , k' q and the length of the elastic L > have to be set such that 
x(L) = X! and /^x(L) = 6\. 



This can be done by means of a shooting algorithm where we use the B-spline solutions (5.441 
(which correspond to the coordinate dependent mode-lines of a product of two Heisenberg approx- 
imations of the Green's functions) as an initial condition. 

The shooting algorithm works as follows: First we write everything in one ODE-system 

/ x ( s ) \ 

I y(s) ^ 
d(s) 
x(s) 

y( s ) 
vi s ) 

k(s) 
\ k(s) 1 
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with initial condition 

(x(0), y(0), 0(0), i(0), y(0), i(0), y(0), k(0), «(0)) 



(0,0,0,1, 0,0,1, Ko ,n' ), 



(7.94) 



where we note by means of left-invariance we can assume that g(0) = e. Now this system of 
equations has a unique solution and can be numerically solved by a standard Runge-Kutta method 
yielding the numeric solution (x(s), y(s), #(s) = arg(x'(s) + iy'(s))). This defines a func tion ip : 



spatial curve s 
on the function <f> 



R 2 ) which maps (kq,k' ) (which determines the initial condition ( 7.94 1 ) to the 
(x(s),y(s)) solution. So finally, we apply a (dampened) Newton-Raphson scheme 



x R x M+ -> R+ given by 
(k ,k' ,L)» ||(^( Ko ,^))( j L)-x 1 || 2 + 7 2 (Z(^( Ko ,^))(L) -e,) 2 , (7.95) 
for suitable choice of 7 > 0, where we use finite difference approximations for the derivatives of 



At this point we note that the initial guess for kq, k' and L, which can be derived from (5.441 
is given by 

K (n) - y"(°) 

(i+(y'(o)) 2 ) 3 /2 • 



«'(0) - 



y'"(Q) 



(l+(y'(0))2)3/2 



■ (y"(Q)) 2 y'(Q) 
(i+(y'(o)) 2 ) 5 / 2: 



L=J ^(y'(x)Y + ldx, 
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Figure 12: A good initial guess, such as (7.961, for the shooting algorithm is crucial as we do not 
want to get stuck in local minima. Typically a dampened Newton-Raphson algorithm k n+ i = 
k n — q £><I>(x„)<I>(x„), < q < 1 in the shooting algorithm, i.e. to find the zeros of (7.951 avoids 
wide jumps, so that the chance of getting stuck in a local minima is reduced. 



(7.96) 

with 2/(0) = , y"(0) = ^(3yi + - 20 o )) and y'"(0) = %(-2 Vl + Xl {6 - ft)), y{x) = 

x &o + ^3- ( — 2t/i + Xi(6o — &i)) + f^(3yi + x\{6\ — 26 ))- To this end we note that this good initial 
guess is highly relevant to avoid the shooting algorithm to get stuck at local minima, as illustrated 
in Figure [12] 

The horizontal curve in SE(2) corresponding to the elastic is given by s i— ► g(s) = (x(s), Zx(s)) 
and indeed satisfies 

5 (0) = (x , e i9 °) e SE(2) and g(L) = ( Xl ,e^) e SE(2). 
Exact derivation of the elastica 

In this subsection we will investigate well-known analytic formulae for (the curvature k(s) of) 
elastica curves, to get some analytic grip on the behavior of these curves. It turns out that 
exact formula for elastica curves involve special functions (Jakobi-clliptic or theta functions) with 
practical disadvantages due to the twice integration of their curvature. 
Consider the ordinary differential system 



(7.97) 



2k" + k 3 = en, 
k(0) = k '(0) = K 'o 

A multiplication of the ODE by k' and integration over arc-length yields 

2^ = -K 4 + 2e K 2 + Ad, (7.98) 
as ) 

where G\ is an integration constant, related to the initial conditions by means of 

From which it directly follows that s is an elliptic integral in n, [43., 
r/, W 2du 

(7.99) 



Ko V-u 4 + 2eu 2 + 4Ci 



which only holds for C\ > (k' ) 2 — V. Now for C\ ^ this can be rewritten as follows 



4 

■»(«) 



Jko VW l -2eu 2 -4:C 1 ■> KO v /(u-7+ ) («+7l ) -£= 
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where v — uj a/t! and where 7^. = e ± Ve 2 + 4Ci are the real zero's of it 1— > u 2 — 2eu — AC\. We 
have 



dv 



where the constant C2 = sign{Cf } i 4 — Jo 



+ • 



dv 



(7.100) 



As a result we can rewrite (7.1001 



k(s) 



71 sn 



^(a + Ca) 7! 
2 '74 



(7.101) 



where sn(-,fc) denotes the Jacobi elliptic function of the first kind, which is the solution of (7.971 
for C\ 7^ if and only if 



d - C\ := ( K > ) 2 + 14 - §k§ and C* 2 = sigu{Cf} / 



%/1-w 2 



7 + 



So we see that the curvature is a periodic function with period 



2 

77+ ./u 



ir/2 



d6 



1 - (1 - ^)sin 2 



7±=c±V c2 + 4C, l £ 

(7.102) 



note that Ci 1— > Tb I)e is a monotonically decreasing differentiable function with lim ^4 — = 0, so 

for applications C\ is typically small since the number of periods over a fixed interval of interest 
corresponds to the number of turns the elastica makes during this interval. For C\ — > solutions 
s 1 ► x(s) tend to straight lines (they are straight lines if ko — = 0). For kq — y/e and k'(0) = 
the solutions are circles. 

Finally we note that the elastica s i— > (x(s),y(s)) = z(s) := x(s) + iy(s) follow by their 
curvature k(s) by means of 



dz 
ds 



; ie(s) and 



d6 
ds 



k(s), 



(7.103) 



now the primitive 8(s) can easily be derived analytically from (7.101 1, but the second integration 
step which provides the actual curve z(s) is a non trivial expansion in elliptic functions, for details 
and derivations see [S3] . This problem is partially resolved in Mumford's approach |15]p.502- 
505, where the Jakobi-elliptic functions are replaced by theta functions and where the arc-length 
parametrization z(s) — x(s) + iy(s) does not involve an integration. But even in this approach 
the standard solutions 

s i-> c^-\og6x(s-ri)-as, with a,c G C and A = Z+itZ.r] = -it/4 or A = Z+(l/2)(ii+l)Z,77 = 0, 
ds 

involve several parameters which are not straightforwardly related (to e > and) the boundary 
conditions x(0) = Xo, x(0) = (cos 6q, sin 9q), x(1) = xj., x(l) = (cos 6\ , sin Q\ ) . So for computation 
purposes a shooting algorithm of the type (7.951 is preferable over an entirely exact approach. 

7.1 The corresponding geodesies 

Next we are going to repeat the proceeding with the "only" difference that we take a square root 
of the integrand so that we have a homogenous energy 

L 



£(x) = J y/K 2 ( S ) + e ds, 
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x G C, which is related to the Cartan connection, recall the second Example in section [6 

We stress that the corresponding lifted curve s i— > (x(s), Zx(s)), recall definition 5.5[ is a 
geodesic on SE(2) in the classical sense, since by our restriction to horizontal curves we have 
k(s) = 9(s) and we can rewrite the energy as 



I V« 2 00+ed S = J y/WsW + mMsWP 2 d*= J \i(s)\ ds = J f(s)^( S )g l3 ds, 

I ,J 



with (3 = e, where we recall (6.701 and (6.781 



The energy after deformation 



7.90 1 becomes 



S(x + hSn)= j \/k 2 new (s) + e ds NE w = J V^ 2 + 2 ^"« + 2<5/ik 3 + e + 0(h 2 )(l - Shnjds 
o o 

L 



o 
o 

L 



hS" K+hSK 3 
K 2 +e 



+ 0(h 2 )^ (l-hS K )ds 

= j 7^+7 fl + hS'^hSK 3 _ §K + 0(h 2 )) ds 

o K c ' 

= £(x) + ft / V>?^ ( feg '^+f 3 - ds + 0{h 2 ) 



where we used a/1 + x = 1 + hx + 0(x 2 ) and (7.91 1. So in order to get a local minima the energy 



should increase under all possible deformations parameterized by i5 and therefore we have 



the solution of which is straightforwardly derived from (7.1041 by substitution 

K 



\/k 2 + e 
which gives us 



(7.104) 



(7.105) 



K 2 ( S ) = 



e(z( S )) 2 



for (z(s)) 2 < 1 



1 - W*)) 2 

with z(s) = ±(z - ^4)e-^ s + §(z + ^)e^ s , i.e. 
z(s) = zq cosh(v / e s) H p sinh(\/e s) 



(7.106) 



with z = / K ° , z[, = — eK ° 3 which is only valid for 



s e [0, s max ) := [0, -= log y - 

V e \ zo + e 2 



(7.107) 



For a comparison between the elastica and the geodesies derived in this section see Fig ure 13 
Here we recall that the e of the elastica curves has to be set to e = (3 2 = 4aDu (recall (5.48 1) 
whereas the e of the geodesies has to be set to 

= PlL 
€ D 22 
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(See (B.152 1 and sec also ( 9.139 1). Both parameters e have the physical dimension [LENGTH] 2 . 

we have set D 2 2 = ~r~ , so that e = 4aZ?i i = . 

zz 4a ' 11 D22 

exact tangible formula (A. 147 1 (where the parameters are given 



In our comparison in Figure 



13 



In appendix [A] we derive an 
by (A. 148)) for the geodesies in the general case by means of symplectic differential geometry and 
Noether's theorem. Moreover in appendix [A] we will derive an important conservation law (the 
so-called co-adjoint orbit condition) along the geodesies and we re-derive (7.104) in a shorter and 
much more structured (but also more abstract) way. 



8 Non-linear adaptive diffusion on orientation scores for 
qualitative improvements of coherence enhancing diffu- 
sion schemes in image processing. 

A scale space representation uj : M. d x R + — > K of an image / : M. d — > K is usually obtained 
by solving an evolution equation on the additive group (R d ,+). The most common evolution 
equation, in image analysis, is the diffusion equation, 

d s u f (x, s) = V x • (C(u f )(x, s) V x u/)(x, s) 
u/(x,0) = /(x), 

where C:L 2 (l 2 xK+)n C 2 (R 2 x R+) -> C^R 2 x R+)) is a function which takes care of adaptive 
conductivity, that is conductivity depending on the local differential structure at (x, s, U/(x, s)). 
In case C = 1 the solution is given by convolution w/(x, s) = (G s * /)(x) with a Gaussian kernel 

G s (x) = — ^-T-e~y with scale, s = \u 2 > 0. 

(4-irs) 2 

As pointed out by Perona and Malik [45], non linear image adaptive anisotropic diffusion (dif- 
fuse less at locations with strong gradients in the image) are straightforwardly taken into account 
for by replacing the isotropic generator A = V • V by V • (c(\\ Vuf(-, s)||)V), i.e. C(w/)(x, s) = 
c(||V x w/(x, s)||), where c : M. + — * R + is some smooth strictly decaying positive function vanishing 
at infinity. This is based on the intuitive idea that if (locally) the gradient is large you do not 
want to diffuse too much. By restricting ourselves to positively valued c > one ensures that 
the diffusion is always forward, and thereby ill-posed backward diffusion is avoided. The most 
common choices are 

c(t)=e W*, c(t) = -_L andc (<) = — 2==, (8.108) 

involving parameters p > 5, c, A > 0. The corresponding flux magnitude functions are given by 

(f>(t) = tc(t), with t = ||Vit/||. 
The sign of 

(j>'(t) = c(t)+tc'(t) (8.109) 

is important, since if <f>'(t) > then the magnitude <j>(t), t = \\ Vw/||, of the flux 

c(\\Vu f \\)Vu f (8.110) 

(by Gauss Theorem) increases as ||VZ7/|| increases, whereas if <f>'(t) < the magnitude 0(||Vu/||) 
of the flux (8.110) decreases as ||Vf7y|| increases. Typically, this introduces an extra "sharpening 
effect" of lines and edges. However, this sharpening effect (besides the decay of the conductivity 
function c : M. + — > M. + ) should not be mistaken for ill-posed backward diffusion because in all cases 
c(t) > for alH > 0. To this end we note that the Perona and Malik equation can be rewritten in 
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Figure 13: Left column; joint plots of the elastica and geodesies in R 2 . Right column; the 
corresponding graphs of the curvature as a function of the y-coordinate of points along the 
curve. We use the same e for both elastics and geodesies. This means we should set e = AaDu = 
j^, where (a,Dn) are the parameters of Mumford's direction process and where (Du,D 2 2) 
are the parameters of the contour enhancement process, also proposed by Citti and SartifTTj 
as a cortical model for contour enhancement. Top row parameter settings: e = 0.0125, Length 
geodesic s — 15, k(0) = —0.20502, k'(0) = —e~i(K 2 ) + e)Ko, so that s max = oo,x = 0, 6> = 0,xi w 
(11.88, -8.43535), 9 1 = -51.8947°. Middle row: Length geodesic s = 15, e = 0.125, n = -0.2, 
k'(0) = -e-i{K% + e)K =*> s max = oo, x = 0,6 = 0,xi w (13.58, -6.09), 6>i = -29.36°. 
Bottom row: Illustration of an extreme case, Length geodesic s — 48.997, k'(0) = —0.0058, k = 
-0.125 «(0) and re'(O) such that x s codcsic (s = 48.997) = x x and Z(x gcodcsic (s = 45), e^) = 6»i, 
x = O,0 O = 0,xi as (—41.383), 6*i = 50.63°. Typically, the geodesies have more curvature at 
the boundaries whereas the elastica curves have more curvature in the middle of the curve. For 
reasonable parameter settings (such as top row where we set e = AaDu = %^ = 0.0125) the 
geodesies are close to the elastica curves. This makes sense since the random walkers in the 
contour enhancement process are allowed to turn in negative ^-direction in contrast to random 
walkers in the contour completion process , recall Figure [3] 
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Gauge-coordinates w along the normalized gradient e w = yvLj y ^ u f ano - v &l° n g the normalized 
vector e„ = nyL II (~ ®v u f> orthogonal to the gradient, using 



.1091: 



ft = divWIlVu/IDVu/) = & (c(fe)fe) +c(fe)? 

_ (h'( du f\ 92 u f J_ c ( du f) d2u f 



(8.111) 



with ^tf = ||v \- jz (V*Uf)H x [uf](V x Uf) T and ^ = ||V X «/||- Now the coefficient in front of 



in the righthand-side of (8.111 1 is negative iff </>'(||Vu/||) < 0, but this does not correspond to 
inverse diffusion. 



The intervals, respective with the choices of c(i) in (8.108 1, where 0' is positive are 

[0,A(2c3>)~£), [0,A(2jj-1) - S*) and M+. 

So the "sharpening effect" does not occur in the case c(t) — , „ . For the other two choices 

there is always a danger that the "sharpening effect" due to switching sign of <fi'(t) can cause 
"staircasing effects", [56] p. 52: That is step-edges will evolve as a staircase over time due to the 
fact that strong gradients will lead to an effective sharpening of the data whereas weak gradients 
will lead to relatively smoothing of the data. 

A further improvement of the Perona and Malik scheme is introduced by Joachim Weickert, 
|57j . who also uses the direction of the gradient V x Uf of Uf, which is not used in the algorithms of 
Perona and Malik type. Therefor he proposed the so-called coherence enhancing diffusion schemes 
(CED-schemes) where the diffusion constant c is replaced by a diffusion matrix: 

S( Uf )(x, S) = (G(j * Vu f (;s)( Vu f (; S )) T )(x) 

C(it/)(x,«) =al+(l-a)e c* I w V )<.,.))-\ a (S(« / >(x,.»)=' e 2 (5(u/)(x, «)) e| , (5(u / )(x, «)) [ ' ' 

where a G (0, 1), c > 0, a > are parameters and where the help-matrix S, with eigen values 
{Xi(S(uf)(x, s))}j=i )2 is used to get a measure for local anisotropy e ( A i( s ("/)<=. s ))- A 2<s(« / )(>=, S ))) :4 
together with an orientation estimate e2(5 , (u/)(x, s)) which is the eigen vector with smallest eigen 
value (orthogonal to the average gradient). In order to get robust/reliable orientation estimates it is 
essential to apply a componentwise smoothing on the so-called "structure-tensor field" Vit/® Vu/. 
The amount of smoothing/averaging of the structure tensor field is determined by a > 0. 

This lead to useful and visually appealing diffusions of the famous Van Gogh paintings and 
fingerprint images, see Figure [14J Nevertheless, this elegant method fails in image analysis appli- 
cations with (almost) crossing lines and contours as it starts to create strong artificial curvatures 
at crossing locations where the gradient is ill-defined. 

As this is a major drawback in many (medical) imaging applications we are going to solve this 
problem by considering similar non-linear adaptive evolution equations on invertible orientation 
scores. To this end we note that in invertible orientation scores crossing lines are nicely torn apart 
in the Euclidean motion group. Moreover, in our orientation scores we have full information on 
both local direction and local curvature (!) at hand which enables us to steer the diffusions in a 
left-invariant manner on the orientation scores (and thereby Euclidean invariant manner on images 
via the unitary wavelet transforms) . See Figure [16] 

8.1 Coherence Enhancing Diffusion on Orientation Scores 

In order to obtain adaptive diffusion on orientation scores we will use the following basic non-linear 
left-invariant evolution equations on SE(2) as a starting point 



/(D11 (U)) (g, t) \ (pd e \ 

d t U(g,t) = ( /3do d S ft, ) (D 22 (U))(g,t) % U ^ 

\ (D 33 (U))(g,t)J \ dj 

for all g G SE(2),t> 0, 
I U(g,t = 0) = W4f](g) for all g G SE(2), 



.113) 
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Linear Non-Linear 
Input image « = i,c = i c(||v x i*/( c\u f ,-,s) 




Figure 14: From left to right: input image / of the well-known portret of Van Gogh, computed on 
comparable slices s) in a linear scale space representation C = 1, Perona en Malik non-linear 
scale space representation (left case in (8.108)) and coherence enhancing diffusion (CED) given by 
([8~TT2| by Weickert, [57]. 



with (3 > (recall 6.701 and where the functions D^k 
C 1 (SE(2) x R+), 



(<M) 



l+), k = 1,2,3 given by 
(D kk (U))(g,t)>0, Ueh 2 (SE(2)x 



h 2 (SE(2) x E+) n C 2 (SE(2) x R+) 



should be chosen dependent on the local Hessian HU(-,t) of U(-,t) (similarly as was done in 
the CED-scheme (8.112 1 ) such that at strong orientations -D33 should be small so that we have 



anisotropic diffusion in the spatial plane along the preferred direction , while at weak directions 
D33 and D22 should be relatively large and isotropic D22 ~ D 33 . Usually we set D22(U)(g 7 t) = 1, 
since in general there is no reason to make D 2 2{U){g , t) dependent on g and in such cases a simple 
re-parametrization of time yields D22 = 1 ■ 
Example 1: 

For example one can take D 2 2(U)(g 7 t) = 1, Du(U)(g,t) = D 33 (U)(g,t) = e~ ' ( ' ' " , where 
c > is a standard (Perona Malik-) parameter where s(J7)(<7, t) is a measure for orientation strength 
like 



s(U)(g,t)=max(-Re(\i(H\U\(g,t))),0), 



.114) 



where Xi(g,t) is the largest eigenvalue of the Hessian H\U(-,t)\(g) = [„4j.4j|£/(-, £)|](g), i = 
l,...3,j = 1,...,3 where i is the row index. Here we stress that we take the Hessian of the 
absolute value |f/(-,i)|, since the absolute value of an orientation score is phase invariant, i.e. it 
does not matter if you are on top of a line (large real part) or on the edge of a line (large imagi- 
nary part), recall Figure [l](d) and recall subsection 6.3.1 In practice we use Gaussian derivatives 
(4.301, rather than usual derivatives, of |t/(-,£)| (so isotropic with scale si in the spatial part 
= and scale S2 in the angular part = T with periodic boundary conditions) at small scales 
si = /3 2 si, S2 > (typically \J2 s\ is in the order of say 2 pixels and \J2 s 2 is in the order of say 
ff ). See Figure 
Example 2: 
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Another modification (or rather slight improvement) in the non-linear diffusion system (8.1131 is 



obtained by replacing the orientation strength s(U)(g,t) (8.1141 in the first example by 



s(g, t) = max{- 



■£(e? (g,t)) T {H p \U{;t)\{9)) T H p \U{;t)\{9) 6?(fl,*) , 0} 

i=l 



.115) 
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-(tt(x + ei, ! + 1) - tt(x + ei, ! - 1) - u(x - et , ! + 1) + u(x - ei.i - 1)) 



-(«(* + 



.' + 1 



I + 1) - u(x + 



,1+1 



Si - i)) 



Figure 15: Finite difference scheme of 



[52) , for sampling on the grid of our moving frame {eg,e^ 
cos9e y }. 



1131 where we use second order B-spline interpolation, 
cos 9 e x + sin 6 e y , e„ = — sin 9 e x + 



where we recall that the symmetric Hessian Hp\U(-,t)\ was given by (6.841 and where e°(g,t) 
and &%{g,t) denote the remaining eigenvectors (with largest two eigen values) of the symmetric 
matrix Hp\U(-,t)\(g)) T ' H \U{-,t)\{g): 



(F /3 |t/(^)|( 5 )) T iJ /3 |C/(-,t)|( 5 ) c(g,t) = X c(g,t) 
(H \U(;t)\(g)) T H[i\U(;t)\(g) e°(g,t) = X 1 e? (g,t) 
(H \U(;t)\(9)) T H \U(;t)\(9) e°(9,t) = \ 2 e°{g,t), 



|A | < |Ai| < |A 2 | 



So the righthand side of (8.115) is to be considered as the Laplacian in the tangent plane orthog- 



onal to the vector c(g,t). Here we recall from section 6.3.1 that c(g,t) corresponds to the best 
exponential curve fit 7(s) 



» E Ci(g)Ai 

ge i=1 to the data |J7(-, t)\. 



With respect to the numerics of (8.113) and (8.1171, we implemented a forward finite difference 



scheme using central differences along the moving frame {9, £, rj} where we used 2nd order _B-splinc 
interpolation, [52], to get the equidistant samples on the {£, ry, 0}-grid from the given samples on 
the {x, y, #}-grid, see figure [l5| thereby our method is second order accurate on SE(2). As our 
algorithm (for details see |28]) is only first order accurate in time, we took small time steps in our 
contour-enhancement experiments. With this respect we note that a Crank-Nickolson scheme for 
time integration is second order in time and can improve computation time. 

Another issue for reduction of computation time is the time dependent conductivity matrix, 
which from a strict point of view needs to be updated at each time step of the evolution. In 
practice, usually the updating of the conductivity matrix in our finite difference scheme does not 
have to be done at every time step and even the linear case where the conductivity matrix is not 
updated at all (so the matrix is determined only by the absolute value of the initial condition 
| [/"(-, 0)| = IWvJI) yields good results. 

8.1.1 Including adaptive curvatures in the diffusion scheme using Gauge-coordinates 



In subsection 16.3.11 we discussed two methods of how to obtain curvature estimates in orientation 
scores. This was done by finding the best exponential curve fit to the absolute value of the 
orientation score (which is phase invariant, recall Figure [l] (d)). We distinguished between two 
approaches. In the first approach we considered the best horizontal exponential curve fit to the 



data (6.85), whereas in the second approach (6.871 we considered the best exponential curve fit to 
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Original 



CED-OS t = 30 



CED t = 30 




Original 



CED-OS t = 2 CED-OS t = 30 CED i = 30 




Figure 16: Medical image applications. Top row 
orientations scores (CED-OS), see ( 8. 1 13[ > and 



1171 



Result of coherence enhancing diffusion on 
and standard coherence enhancing diffu- 



sion directly on the image, see ( |8.f 12 ), (CED) of bone-tissue. Middle row: Result of coherence 



enhancing diffusion on orientations scores (CED-OS) and standard coherence enhancing diffusion 
directly on the image (CED) of 2-photon microscopy images of a muscle cell. Bottom row coher- 
ence enhancing diffusion on orientation scores(CED-OS) and standard coherence diffusion (CED) 
on medical images of collagccn fibers of the heart. All these applications clearly show that coher- 
ence enhancing diffusion on orientation scores (CEDOS) properly enhances crossing fibers whereas 
(CED) fails at crossings: CED creates a "van Gogh" type of painting out of the original image /. 



the absolute value of the orientation score. Both approaches yield a curvature estimate which in 
this paragraph we assume to be given. We shall write (n est (\U\))(g,t) for the curvature estimate 
of the score U via its absolute value \U\ at location g £ SE(2) at time t > 0. Since we only want 
to include curvature at strongly oriented structures we shall multiply it with a front factor: 



K(\U\)(g,t)= (l-e V-33( S . t) ; j Kest (\U\)(g,t), 
where d K controls the soft threshold on including the curvature estimate K est (\U\)(g,t). 



Now we can include curvature in our scheme (8.1131 by replacing i— ► + ndg. To thi s end 



we recall that the exponential curve s \— > e s ( 9 e+ K °e'U = e sd x+Kd e yields a circular spiral (6.731 



whose projection on R 2 is a circle with radius \k\ 1 if n is constant. Moreover, along horizontal 
curves we have 

-C/( 7 ( S )) = k(s) — (~/(s)) + -gr(l(s)) 



where n(s) — || and (d£, 7(s)) = 1, see Appendix [c| See Figure 
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Here we should be careful 
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Figure 17: Illustrations of the heat-kernels : SE(2) ->■ R+ on SE{2). Left: D = 

diag{_Dn, _D 2 2, 0} in left-invariant coordinate frame {dg^d^d^}. Right: D = diag{0, Dbb, 0} in 
gauge-coordinate frame {d a ,db,d c } with 7 = 0, n = 0.06 ,(3=1, Dbb — 1.0036 and t = 70. 




Figure 18: The gauge coordinate frame {d a ,db,d c } illustrated with respect to the basis of left- 
invariant vector fields {dg,d^,d v }. Here we note that the curvature estimation k is given by 
(6.871 and {d ai db,d c } are g iven by ( |8. 1 16 ) , where db is the direction determined by the optimal 
exponential curve fit |6.8l| to the data The angle 7 = j(U)(g,t) = axg{d(g,t) + icl{g,t)) 
intuitively tells us how "horizontal" the orientation score is at location g £ SE(2) at time t > 0. 



since {dg,d^ + kc^c^} are (in contrast to ^{/3dg, <9j, d^}) no longer orthonormal with respect to 
the (•, -)/3 inner product. Therefore we are going to introduce the gauge coordinates, aligned with 
the optimally fitting exponential curve 

3 

s^g expCs^cKs,*)^), c*(g,t) = (cl(g, t),c\{g, t),c^(g, t)) £ R 3 , 

i=l 

with Hc*!!^ = (c®) 2 + [3 2 {c\) 2 + /? 2 (c*) 2 = 1, to the orientation score data \U(-,t)\ at position 
g £ SE(2) at time t > 0. These Gauge coordinates are given by 

< d b = /3{cid s + gdr, + cld e ) 

d c = 1 1 c * =df H — , 5 c * = G>n . 

Note that the gauge-vector is along the best exponential curve-fit direction, i.e. db = c* and note 
that the span of the tangent vectors {d a , db} corresponds with span{<9 Q , d c } = (c*)- 1 . For geometric 
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Figure 19: Illustration of the projection P^eb of the vector field = db on the spatial plane 
(image-plane) plotted on fixed orientation layers |E/|(-, 9) of the absolute value \U\ of the orientation 
score U. Along this vector field we plotted circular arcs to also include the curvature k. These 
circular arcs correspond to the projections of best exponential curve fits to the data \U\ at each 
location g on the image plane. From left to right, the original image, and plots for 9 = ^ and 
9 = — ^. Note that at positions g € SE(2) near strongly oriented structures, the vector is 
better aligned with the local image structure than the vector. 



understanding it helps to consider the Gauge tangent-vectors in ball-coordinates with respect to 
the basis of left-invariant vector fields {dg,d^,d v } so that it becomes obvious which rotation in 
50(3) (or rather which class of rotations in 50(3) / 50(2) = 5 2 , if we do not distinguish between 
directions in plane (c*)^) is required to map the standard left invariant basis {dg, d^,d v } into the 
basis Gauge-coordinates. See Figure [T8| The Gauge-coordinates in ball-coordinates read 

d a = — cos a cos 7 d^ — cos a sin 7 d n + sin a dg , 

db = sin a cos 7 d^ + sin a sin 7 d v + (5 cos adg , (8.116) 
d c = — sin 7 d^ + cos 7 d v , 

where the Euler-angles read 

a = arccos cf = arccos , K , 

7 = arg(cf + i c'J). 

Here the function 7 which maps U to j(U)(g, t) = arg(cf (<?, t) + i c*(g, t)) See Figure 19 

Now the diffusion generator diagonal along the left-invariant Gauge vector fields is given by 

D aa (d a ) 2 + Db b (d b ) 2 + D cc (d c ) 2 = ( Pd e d i d ri )Ml 



D aa 





x 




( (3d 



























I sin a — cos a cos 7 —cos asm 7 \ 

where jf^, = cos a cos 7 sin a sin a sin 7 is the rotation matrix in 50(3) which maps the 

y — sin 7 cos 7 I 

tangent vector (3dg to the tangent vector db- 

Now by straightforward computation this leads to the following non-linear evolution equations 
on orientation scores 

' d t U(g,t) = ( fide d t d v )^^x 

D bb K 2 + D aa f3 2 KtJ(O bb -I) aa )cas-y K0 <.r> bb - D aa ) rfn 

««%-D IO )™7 O cc (k 2 +/3 2 ) + {(D bb -O cc )/3 2 + (D aa - D cc ) K 2 ) cos 2 7 cos 7 sin 7 ( K 2 (D aa - D cc ) + ft 2 (D bb - D aa ) ) 

K 0(D bb -D aa ) sin 7 cos 7 sin 7 ( K 2 (D aa -D cc ) + I3 2 ( D bb - D cc ) ) D bh f3 2 + D aa K 2 + ((D cc - D bb )p 2 + (D cc ~ D aa ) K 2 ) c 

/ /39o \ 

x 3« \U(g,t), for all 9 £ SE(2), t > 0, 

V a, / 

E/(ff,t = 0) = W^[f]{g) for all g 6 S£(2), 
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3.117) 



where for the sake of clarity we used short notation Da — (Du(U))(g,t), for i = a,b,c. Now 
again we set 

D bb = 1 and (D aa (U))(g,t) = (D cc (U))(g,t) = e~ """l"'*"' , c > 0. 



Here we take (8.115) as a measure for orientation strength. In the Gauge coordinates this measure 



can be written 

s(g,t) = max{-A c x|£/(-,t)|( 5 ),0} = max{- ((d a ) 2 \U(; t)\ + (d c ) 2 \U(-,t)\) (g),0}. 



and the conductivity matrix in (8.117) simplifies to 



ii 2 



/ k 2 + Daa/3 2 k/3(1 -D aa ) cos7 K0(l-D aa ) sin7 

— K /3(l-D aa )cos7 D aa (n 2 +l3 2 ) + (l-D aa )/3 2 cos 2 7 cos7sm 7 /3 2 (l-D aa ) 

1 M 1 COS 7 sin 7 /3 2 (l-D aa ) /3 2 +£>aaK 2 + (Daa -l)/? 2 COS 2 7 



Kf3(l-D aa ) sin 7 



See Figure [17] for an illustration of the special case D bb is constant, Daa — D cc = 0, 7 = 0, which 
despite the strong degree of degeneracy still leads to a smooth and useful Green's function since 
the Hormander condition, recall subsection |4.2.1[ is satisfied. 

Remark: 

Although not discussed here it is worthwhile to consider the components {C v (U)(g, s)} of the 
adaptive conductivity matrix with respect to the basis of left-invariant vector fields {Ai, A2 , A3} := 
{de, d^,d v } on SE(2) as the inverse-components of a metric attached to the graph of g 1— > U(g, s). 
So rather than changing the constant conductivity by an orientation score adaptive conductivity 
within the di ffusion equation (as we did in this section) we could change the constant left-invariant 

3 3 

metric (6.70) by an orientation score adaptive metric G(U)(g,t) — Gij(U)(g,t)dA l ® dAP, 

i=ij=i 



where Gjj{ U){g, t) are the components of the inverse matrix of [G 13 (U) (g, t)] := [G 13 (JJ)(g, t)]. So 
in stead of (8.117) we could consider the Laplace-Beltrami flow 



E EA{ y/det G(U) (g, t) G ij (U) (g, t) Aj U(g, t) } (g, t) for all g 6 SE(2), t > 0, 
(=1 j=i 1 > 



d t U(g,t)= , 1 

{ U(g,t= 0) = WAf](g) for all g G SE(2). 

This only means (by the product rule for differentiation) that we must add to the righthand side 
in the PDE in (|8.117h which equals 



EE- 4 ' i Gij ( U )(9, t)AjU(9, t)} (g, t), 

i=l 3 = 1 



the following terms: 



Vdet G(U)(g,t) (z 



J2 Ai\yJdetG(U)(g,t) \ (g,t) E G i3 (U) (g, t)Aj U (g, t) . 
i=i *• > j=i 



9 Towards a graphical eraser: Morphology PDE's on SE(2). 



Sofar we considered automatic line and contour enhancement via linear and non-linear left- 
invariant diffusion equations on (invertible) orientation scores of 2D-images. In this way we ob- 
tained an automated "graphical sketcher" . What is missing is an automated "graphical eraser" 
which erases brush strokes which are to far away from the the zero crossings dgU and d v ll and 
which make the completion fields to broad. See Figure [20j 
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Figure 20: The sections [T][2j [3] |4][6j|8] provide the theory for automatic sketching process by means 
of left-invariant evolution equations on invcrtiblc orientation scores, where the general scheme is 
explained in Figure [4] The sections [5] [7] [5] and Appendix [A] [S] [C] serve the practical purpose 
of developing an automated erasor by means of left-invariant curve extraction. Section [9] serves 
the practical goal to create a continuous process from a sketched image to an image consisting 
only of these curves, where we again use the general scheme explained in Figure [4] but where 
we replace the (non-linear) left-invariant evolution equations (8.1131 and ( 8.117[ > (in sec tion [8 1 by 
Hamilton- Jakobi equations (also known as morphology equations in image processing) (|9. 118 1. 



Moreover, our oriented wavelets are not "perfect" line detectors. They often yield too broad 
distributions 9 i— ► f//(x, e l6 ) at fixed positions x G M 2 where a line/contour is present. Therefore 
it is often desirable to erase (i.e. to apply erosion) in both 9 and rj direction. 

In the previous chapter we have seen that the modes of the direction process are for reasonable 
parameter settings (^n small) closely approximated by the intersection of the smooth surfaces 
{g G SE(2) | d e U(g) = 0} and {g G SE{2) \ d n U(g) = 0}. The gradient V 2 U = (d g U,d v U) of a 
(diffusecp^) orientation score U G L,2(SE(2)) locally points to these zero-crossings. 

Therefore we propose the following 2 PDE-systems on real-valued (processed) orientation 
scores U : SE(2) -> M 

d t W(g,t) = -\\{CV 2 )W{g,t)\\ 2Q = -((d) 2 (WM 5 , i)) 2 + (csW^,*)) 8 )* 
W(g,0) = \U(g)\, with a = D u and c 3 = D 33 , 

C G with C = diag(ci,C3). For £ — 1 the solution is given by 

W(g,t) = (he SE (2) U)(g) := inf [U(h) - k^{g- l h)\. (9.119) 

n,£i>g 

where the morphology kernel k t is given by 
K DllDaa (n) 

k t {g) = * n {y> (9.120) 

where K® 11 ' 033 (g) equals (after rotation by ir/2 in each spatial plane, mapping £ to 77) the diffusion 
kernel on SE(2) studied in section [4] Theorem |4.1| (with analytic approximations in section 4.2 1. 
Here we note that the corresponding dilation operator is defined by 

(K® sm U){g):= sup [K^g) + U(h)] 

h£SE(2) 



which is the equivalent of a 5i?(2)-convolution, (2.11 ), where we replaced the (+, -)-algebra by the 
(max, +) algebra. 

19 Sincc thinning and erosion algorithms can be useful both before and after a diffusion step, we simply write 
U 6 I.2 (SE (2)). Here U could be U = |W^/| or U = <£>t([W,/,/|), where t is the stopping time of the non-linear 
left-invariant diffusion or U could be the completion field U = (A — "fl)~ 1 \W^,f\(A* — 7-f) - 1 |W^,/| obtained by 
linear left-invariant diffusion. 
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Remark 9.8. Due to the non-commutative nature of SE(2) (and Hormander's theorem, [36) 
both the diffusion kernel in Theorem 4.1 and the corresponding kernel k t given by (9.1201, in 
contrast to their well-known analogues on IR 2 : 



K t (x,y) = ±- t e-^ and k t (x,y) = log^^ = °^-,(x,y) €R 2 , (9.121) 

are not separable along the exponential curves, respectively, in the (+, -)-algebra and (max, +) 
algebra. For if a differentiable kernel K is separable K(a, b) = k 1 (a)k 2 (b) (resp. k(a, b) = fci(a) + 
k2(b)) and isotropic (adb — bd a )K(a, b) = (resp. (adb — bd a )K(a, b) = 0) then clearly it must be 

a 2 +f> 2 

equal to K(a, b) = /ie * and k(a, b) — 7 (a 2 + b 2 ), for some separation constants /i, A, 7 > 0. 

Before we motivate our conjecture we first explain why the standard approach on finding 
viscosity solutions of morphology equations on M 2 by means of the Cramer transform is not 
easily generalized to SE(2). The homomorphism between dilation/erosion and diffusion/inverse 
diffusion is given by the Cramer transform C = 5° log °£, [2], |10j . which is a concatenation of the 
multi-variate Laplace-transform, logarithm and the Fenchel transform (mapping a convex function 
c : R" — * K onto the function x 1— > [3c] (x) = sup y [y • x — c(x)]). This is due to the fact that 

C(f * g) = $log C{f * g) = £(log Cf + log Cg) = 3 ^log Cf 3 log Cg = Cf® Cg. 

Now, since SE(2) is non-commutative the irreducible representations are no-longer 1 dimensional 
and the Fourier/Laplace transform on SE(2) becomes relatively complicated (sec [19 )App. B). 
Moreover, it is not clear how the Fenchel transform should be generalized to SE(2). 



9.1 Morphology and Hamilton- Jakobi theory 

In the calculus of variations and the corresponding Hamilton-Jakobi theory on some finite di- 
mensional manifold G, with local coordinates {<?i}f =1 , one usually starts with a Lagrangian 
L : M. + x G x T(G) — > M. + of class C 2 , giving rise to the following energy on curves 



£(j)= J L(t )7 (t),Y(f))dt, (9.122) 

to 

where, for now, we shall assume non-degeneracy of L with respect to its dependence on T(G), i.e. 



det (V 2 L(t,g,g)) ^ for all t > 0,g G G. (9.123) 

so that we can express g l as functions of the canonical variables (t,g 3 ,Pi). 

^ = 0(4,^', Pi), (9.124) 

where the components momentum p = V g(L(t, g, g)) are given by pi — dLtyt Qgi ^ ■ 

Here one should make a clear distinction between the parameter dependent non-homogeneous 
case (recall the elastics in section [5] where traveling time coincides with spatial arc-length, t = s 
and where G = SE(2), L(s,[i\,tf]) = ((9(s)) 2 + e)\\x(s)\\ = (6(s)) 2 + e) and the parameter 
independent homogeneous case, recall the geodesies in section [7T] with 



G = SE(2)/Y and L(s, [ 7 ], [7]) = yf {0(s))* + e\\x( S )\\* . (9.125) 

Moreover, one should make a clear distinction between the cases where t > is an intrinsic 
coordinate on the manifold and the cases where time t > is an independent coordinate, see 
@5]p.44-p48. 
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For now we will consider the non-homogeneous case with time as an independent coordinate. 
Soon we will consider t > as an intrinsic coordinate on the manifold, which will turn out to 
be important to relate the direction process to morphological equations on SE(2), whereas the 
relation to contour enhancement processes and morphology on SE(2) requires t > to play an 
independent role. 

Theorem 9.9. A 1 -parameter family of 'hyper -surfaces S(t,g) — a is geodesically equidistant with 
respect to a non-degenerate Lagrangian, that is dot (\7 2 g L{t, g, g)) ^ and L(t, 7 (t), j'(t)) — j§r(i) 
for all congruency curves 7 (satisfying V g L{t, ^(t), "f'(t)) = V g S(t,j(t))) if and only if 

dS 

~ + H(t,g,VS) = 0, (9.126) 

where the uniquely corresponding non-vanishing Hamiltonian function H G G 2 (R + xGxT(G)', M. + ) 
is given by 

H(t,g,p) = -L(t,g,<f>(t,g,p)) + (p,<f>(t, g,p)). (9.127) 
In particular the characteristic function W : R + x G given by 

W(t, g; t Q , go) = inf{ / L(t, 7 (t), 7 (t)) dt \ 7 G C 2 (R + , G) with 7 (t) = g, j(t ) = go}, go e G, t > 0, 



satisfies the Hamilton- Jakobi equation. Along congruency curves the following fundamental equa- 
tions hold 

p = -V g H(g,p) and g = V p H{g,p). (9.128) 

For proof and more background see [l6]p. 12-25. 

Theorem 9.10. The Hyper-surface 

{g G G I S(t,g) = ct + R}, R>0, (9.129) 

where S(t, g) is a solution of the Hamilton- Jakobi equation, is the envelope of the set of geodesic 
spheres 

S (go .t o) , R {^9) G M + x G I W{t,g;to,go) = R} (9.130) 

of radius R > centered on the hyper surface {g G G | S(t,g) = a}. If moreover, the Weierstrass 
excess function E : R+ x G x T(G) x T(G) -> R given by 

EU.g '.<,'.;,'< = L{t,g\f) - /.('.//'. /r) - (g 3 - 9 3 ) dL % g -; g]) > 

is strictly positive, where if denotes the tangent vector an arbitrary extremal curve T issuing from 



go to an arbitrary point g G G such that 
f 9 L(t,g j ,-g J )dt = R, R>0, 



we have that the hyper surfaces ( 9. 129\ ) are supporting hyper surfaces of the geodesic spheres 
centered on the hyper surface {g G G | S(t,g) = a}. 



This results puts an analogy between morphological PDE's, which are of the type (9.127), and 
Huygens' principle. Now consider G = SE(2) and put the following Cartan connection on SE(2) 
given by uj(X g ) — d£(X g )d x , X g G T g (SE(2)), this means that the horizontal part of the tangent 
space is given by the span {dg, d v }. Set C(s, [7], [7']) = j^(^) 2 + iB^(s?) 2 and we § et 

H( P i, P3 ) = -L(2D llPl ,2D 33P2 ) + 2D llP i + 2D 33 p\ = D llV \ + D 33 p 2 2 
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and thereby the morphological PDE-system (9.131) coincides with the Hamilton- J akobi equation 



(9.1271, with £ = 1, and by the above theorems (using the moving frame {S,,rj,9} of reference to 



ensure left-invariance) we see that the viscosity solution is indeed of the type (9.1191. Similarly, if 
we consider G = SE(2), with Cartan connection uj(X g ) = df](X g )d y , X g G T g (SE(2)), this means 
that the horizontal part of the tangent space is given by the span {dg, d{\. 

H(pup 3 ) = -L(2D llPl ,2D 22 p 2 ) + 2D n pl + 2D 22 p\ = D n pl + D 22 pj 

and Hamilton- Jakobi system 

8tW(g,t) = -D 11 (W e (g,t)) 2 - D 22 (W e ( 9) t)) 2 , 
W(g,0) = U(g), ,geG,t>0 



It is still the question though whether the structure element (i.e. the max-plus convolution kernel) 
h{g) = jk(d), with k(g) = W(l, g,0, g ), indeed satisfies k t (g) = -log K ' Dll , fl2i |^ ? 

Considering the usual K 2 -case, i.e. morphological systems on images, this relation does hold 



( |9.121| and it is the question whether the commutative nature of M 2 plays a crucial role. The 
next results indicate that the non-commutative nature of SE(2) does not cause problems in this 
respect. In fact they subscribe our conjecture. 

9.1.1 Hamilton- Jakobi theory and (lifted) elastica curves 

Recall from section [5] that the elastica curves correspond to the modes of the direction process. As 
in a direction process a unit speed grey-value particle is always suppose to move in ^-direction it 
makes sense to identify the arc-length s > with £. By doing this we use the temporal parameter 
t = s = £ as a spatial parameter in SE{2). As pointed out in [3B]p. 44-48 this requires a different 
approach in Hamilton Jakobi-theory. 
We rewrite the Lagrangian 

L(s, B{aU(8),r,(a),e{aU{a),f,{a)) = T (-^((?( s )) 2 + a- 1) = ((9(s)) 2 + e) (d& x(a)), (9.132) 



with (d£,±(s)) = ||x(s)|| = 1, in (9.1221 as follows 



L*(\g(a)], [</(*)]) - L*(6(a),a, V (a);9(8),i(a),f,(a)) ■ £(s) := 
L(8,6(a),T,(»),i$,$)'fa), 

and therefore the components of the canonical impuls p* of the new Lagrangian L* equal 

* _ dh* _ dL _ -j- e _ 

Pl ~ ae 7 ao ~ + 2Dn - Pi' 

pI = l - e(s)pi 

P*3=P3 = 

and the Hamilton- Jacobi equation on G = SE{2), where we again restrict ourselves to horizontal 



curves and horizontal subspaces spanned by {dg,d^}, corresponding to the Lagrangian (9.1391 is 



H*{[g\,p) =p 2 + H([g],p) =p 2T (D n pl - a) = (9.133) 
So the corresponding Jakobi-Hamilton equation is given by 



(9.134) 
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If we now drop our identification between £ and s this yields the following non-linear morphology 
process on orientation scores: 



aw dw i n /w\ 2 
"S3" - "5T ± ^11 I "ST J 

W(-,s =0) = U f 



(9.135) 



which is the morphological equivalent of the forward Kolmogorov equation (2.12) of Mumford's 
direction process (2.131. 

For the corresponding Heisenberg approximation of Mumford's direction process, where the 
curve-length s parameterized is replaced by x rather than £ and where the Lagrangian is simply 
given by L(x,y(x),9(x)) = (9'(x)) 2 + e, 9(x) = y'(x), one can follow the same scheme. Following 
46jp. 44-48 we for the moment introduce an independent time variable r, with x'(t) ^ 0. Later 
on we set r = x 

L(x,y(x),6(x)) = {9{x)f + e, 9{x) = y(x), 

L*(x(r),y(r),9(r)) = x'(r) L(x(r),y(r),9(r), 

canonical impuls vectors (i = 1 : 6, i = 2 : x, i = 3 : y), recall e = AotDn 



* _ dL^ _ 

Pi — d6' — 

* _ 8L^ — 
P2 — dx> — 



i = 



L-piy(x) -p 2 9(x) 



and consequently the Hamiltonian is given by 



H(x,y,9,px,p 2 ) = -L + y(x)p 3 



? 2 

Pi _ Pi 



3 - e 



2 4 

yielding Hamilton- J akobi equation: p 2 + H(x,pi, p 3 ) — 



1 



1 



d x S + 9d y S = + 7 (d e Sy + e<^A 2 S = -{A^Y + e 



(9.136) 



which is again rel ated t o the exact case (9.133 1 by replacing cos 9 by 1 and sin 9 by 9, i.e. replacing 
Ai by Ai, recall (4.361. The characteristic function (computed by the B-spline solutions, recall 
(5.441 now equals : 



I((v"(r)) 2 + e)dr \ y(0) = 0, y'(0) = 0, y'(x) = 9, y(x) = y 
o 



ex = Dn log 



Kf' Dll (x,y,9) 
l (x,0,0) 



A". 



The canonical equations ( 9.128[ ) for the congruency curves through the geodesically equidistant 
surfaces {(x,y,e te ) £ SE{2) \ S(x,y,9) = <r} in the Heisenberg-approximation are now given by 



dpi dH*(x,p) 

~dt ~ dx~ ' 

dxj dH*(x,p) 

dt dp* 



x'{x) = 1 and 

y'{x) = 9{x) 



Pi{ x ) =Pa(x) 
p' 2 (x) = 
p' s (x) = 



with xi = 9,x 2 = x,x 3 = y, from which we directly deduce — » y""{x) — Q,9{x) — y'{x) yielding 
indeed the B-spline solutions (5.441. The impuls vector p along a B-splinc mode starting at 
(x , y ,9 ) ending at (x 1 ,y 1 ,-9 1 ), x < x\ is given by 

' 6yi+2xi(fli-2e ) _ 6x(2y 1 +x 1 (e 1 -e )) 



p{x) 



1 

2Di 



Pi{x) 
Pi{x) = 1 

„ („\ - 6(2 i/1 +x 1 (8 1 -8o)) 

P3\ X > - 2D^xJ 
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In the Heisenberg approximation case we can compute the characteristic function by means of the 



S-spline modes (5.441 yielding 



W(x, y, 9; 0, 0, 0) = min / ^-{ y "{ x )fdx \ y(0) = y, y(0) = 0, y'(Q) = 0, y'(0) 



.0 

_ 3y 2 +3xy6+x 2 8' 2 
~ Dux 3 



which can easily be checked to satisfy the Hamiltonian equation (9.1361. 

Finally, we note that the structure element of the morphology related to the direction process 
is again related to the direction process resolvent Green's function by means of 

W(x, y, 9- 0, 0, 0) = -4 • log ( ffZ^'nT l ■ ( 9 - 137 ) 



Remark 9.11. Consider the corresponding PDE-systems for a direction process and morphology 
on W 2 , then we again have 

{d x -Ddl-aI)R a , D {x,y)=a5 0fi -> R a , D (x,y) = -^ e -^--fe x (y) = -log = £, 

and k x (y) = ^ is the viscosity solution of the Hamilton-Jakobi equation d x k(x,y) = dyk(x,y), 
with fc(0, y) = Sq, with Sq the denotes the convex Dirac function given by <5§(y) = oo if y ^ and 
5 c (y)=0i£y = x. 

Before we proceed with the homogeneous case we put a relevant observation on the exact case. 
Remark 9.12. Despite the fact that k(s) — 6(s) for horizontal curves, the exact elastica equation 



(7.921 is not of the standard Lagrangian type d 9i L + d t {d g 'C — d t d g "C} — 0.. The reason for 
tins is that the arclength parameter s > is curve dependent. Recall that in the only relevant 
horizontal pertubations are: 

X-new(s) = x(s) + hSn(s) 
0new{s) = 0(s) + arctan ^ffv 

we have dSN d ^ w (s) = 1 — h5(s)9(s) + 0(h 2 ). Therefore in order to embed elastica curves in 
the standard Euler-Lagrange/Hamiltonian theory we must follow [53]App. C and rewrite the 
Lagrangian in a curve dependent way 

where s'(t) = ||x'(i)|| = yfWWF+JrfW^ and 

in our notation we distinguish between 9'(t) = 
j- t 0{s{t)) and 6(s) = £0(s). The curvature along the curve t i-> g(t) = (f(t) cos(0(f)) - 
jj(t) sin(0(i)), £(t) sin(0(t)) + r)(t) cos(0(t)), 0(i)) at time t > equals 

K(s(t)) C(tW'(t)-C'(tW(t) 

and consequently 

d e s' = £ and d n ,s' = %>s' = d, r s' = 

%K = -2^2^' + jfrjsrf, d^K = - jfjyCi ^»K = -j^3, d v »K= ~(y)3- 

and thereby as pointed out in |53]App. C the standard Euler Lagrange equations yield 
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Now for horizontal curves Z±(s(t)) = Zx'(t) = 6(s(t)) one can also express the curvature along 
the curve at time t > 



«(s(t)) = 6(s(t)) 



s'(t) 



yielding the equivalent angular Euler-Lagrange equation 



ediL - 9L 



els'- 



{dgL} = 2k(s) + k 3 (s) + en(s) = 0, with k(s) = 6(s). 



9.1.2 The homogenous case: Hamilton-Jakobi theory and geodesies in SE(2) 

Sofar we restricted ourselves to variational problems based on a non-degenerate Lagrangian, i.e. 
a Lagrangian satisfying (9.123) on SE(2). For the elastica curves studied in section [5] this was 
fine, however for the geodesies studied in section [TT] it is not. With this respect we note that if 
a Lagrangian L{g,g) is homogeneous, that is if L(g,Xg) = \L{g,g), then "VgL(g,g) ■ g = L(g,g) 
and thereby det(VgVgL) = 0. So we can not express g in the canonical variables g,p like we did 
in (pl24|! 



Therefore the non-degeneracy condition (9.123) is replaced by another non-degeneracy condi 



tion: 

det(gij(g,g)) ^ 0, 

where the fundamental tensor gijdg 1 ® dg 3 is given by 

l d 2 L(g : g) 
9ij{9, 9) - 2 a p Q g ■ 

and the canonical variables p = Pidg 1 in the non-homogeneous case are now replaced by the 
following variables 

V = 2AoV , where y. L = (g,g) g j = L(g,g) 9Ll ' ' ' ' 



dg 1 



allowing us to rewrite the energy in 



£ (7) = J y / .9ij(7(s),7(s)) f{s)i j {s) ds, 

so 

where gij(g, g)gjk(g, 9) = ^fci which is the length integral in the manifold G in Riemannian 
imiltonian now reads 

(9.138) 



geometry. The Hamiltonian now reads 

H 2 (g k ,y h )=g ij (g k ,yh)yiyj 

and consequently we get 



9 \9 ) (Jh) — 2 dyidyj ' 



g* = H{g\ 9ij (g, g) ^^l^tMl = //,;.,;-. Vl) ""J 

yi=9ij{g,g) g 3 

from which we deduce the following relation between Hamiltonian and Lagrangian: 

H(g h , Vh) = L(g h , g h ), with y h = g hj (g, g) g j . 
As a result, for details see |46jp. 166-170 the corresponding Hamilton Jabobi equations now become 

H(g,VS) = ±l. 



9H{g\ yi ) 



GO 



Now we return to our case of interest, where we equip the principal fiber bundle Px = (SE(2), SE(2)/X, it, R), 
X = {(x,0,0) i € 1} and where we recall n(g) = [g], R g h = hg set connection oj(X g ) = 
(d£, X g )d x so that with this convention the horizontal part of the tangent spaces is given by 

spaa{fy,$,} CT(SE(2)) 

on which we apply the homogeneous Lagrangian: 



L(s, 9(s),r,(s),9(s),f,(s)) = + ^(^)) 2 , (9.139) 

whose corresponding geodesies were studied in section |7.1| The corresponding Hamiltonian is 
given by H{[g\,y) = \J -Dii(yi) 2 + Diiiyi) 2 and the Hamilton- Jakobi equation now reads 



i=Wf^V a5X2 



89 ) \dr) i 

now by setting W(g,t) = S(g)P(T = f ) = ae~ at S(g) we get the following morphological system 

W=^u(^) 2 + J D 22 (^) 2 
W(g,0) = U(g) 



which is exactly (9.118) for ( = \ 



2 ' 



9.2 Graphical thinning 

Another approach for narrowing down the orientation scores (and/or the corresponding completion 
fields) around the zero crossings of dgU and d v ll of some real-valued function U : SE(2) — * K 
is what we call "graphical thinning" . Here oriented grey- value particles are transported to the 
modes, rather than being erased. We propose the following two models for graphical thinning: 

dtW(g,t) = -CV 2 • (W(g,t) CV 2 W{g,t)) = -c 2 (W(g, t) W ee (g, t) + W$(g, t)) - (W(g,t)W vv (g,t) + W%(g,t)), 
W(g,0) = U(g) 

(9.140) 

and its linear counterpart 

f d t W(g, t) = -CV 2 • (W(g, t) CV 2 U\)(g) 

{ = -c 2 (W(g,t) U e o(g) + W e (g,t)U e (g) - (W(g,t) U vv (g) + W v (g,t)U v (g)), (9.141) 
{ W(g,0) = U(g). 

These PDE's have the advantage that the total Li-norm is preserved. However, these equations 
are unstable near the zero crossings of dgW and d v W, which causes serious problems in practice. 
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A Derivation of the geodesies by means of reduction of Pfaf- 
fian systems using Noether's Theorem 

Next we apply Bryant Griffiths approach [9] on the Marsden-Weinstein reduction for Hamiltonian 
systems |40j admitting a Lie group of symmetries on Euler-Lagrange equations associated to the 
functional f v /k 2 (s) + eds, to explicitly derive the solution curves s h- > 7(5) in SE(2). Recall that 
in section 



7.1 



we derived the curvature of the minimizer of J \J n 2 (s) + eds by solving an ODE 
for k that we derived from Euler-Lagrange minimization. Here we will derive the same equation, 
taking into account the restriction to horizontal curves, in a much more structured way, avoiding 
extensive computations, by means of symplectic geometry. Moreover we will derive an important 
underlying conservation law and by the Marsden-Weinstein reduction we will derive the curves 
themselves. Similar to section [7] we will not restrict ourselves to curves of fixed length. 

Consider the manifold Q — SE(2) x M+ x M x R with coordinates (x, y, e l9 , a, k, t), where 
a = ||x'(t)|| so that ds = adt. On Q we consider the Pfaffian equations 



= d£ — adt = 0, a > 0, £ = x cos + y sin 0, -q = —x sin + y cos 8, 

= d?7 = 0, (A.142) 

= d0 - nadt = 0, 



note that these Pfaffian equations uniquely determine the horizontal part I(Q) of the dual tangent 
space T* (Q) , where we recall that along horizontal curves we have 4^ = c _1 4f = K, (dr/, x'(t)) = 0, 

(d£,x'(t))=<7. _ " 

We would like to minimize the energy J \J k 2 + eadt under the side conditions (A.142 1, then 
the gradient of the energy should be linearly dependent on the gradient of the side condition and 
therefor we set 

tp = \/ k 2 + eadt + Ai(d0 - Kcrdt) + A 2 (d^ - adt) + \ 3 dr] 

where Ax, A2, A3 are Lagrange multipliers. Formally speaking, we consider the affine sub-bundle 
Z = { Z q \ q e Q} = Q x T(SE(2))* of T*(Q) determined by 



Z q = { VK^eadt] El q C T*(Q)}, 



3 



Z = Q x T(SE(2))* by the isomorphism (q, A) <-> \/k 2 + eadt\ + E A 



At 



afcl 



fe=l 



Next we compute the exterior derivative of tp 



dip = Vk 2 + ed0 A d< + ^=dK A dt + A 2 d0 Adi; + dA 2 Ad{- dA 2 A adt - A 3 d0 A d£ 
— A 2 dcr A dt + dA 3 AdT] + dAi A d0 — nadXi A dt — erAidre A dt — KAidcr A dt 

where we used the following two equalities 

dd£ = d(cos 0da; + sin 6dy) = - sin 0d0 A da; + cos 0d0 A dy = d6 A (- sin 6dx + cos Ody) = dO A d?7, 
dd^ = -d0 A d£. 

The exterior derivative dip determines the characteristic curves (in our case the geodesies) by 
means of 

i(t)\dip l(t) = 0, and 7 *dt ^ 0. 
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(A. 143) 



So the Pfaffian equations for decent parameterizations satisfying 7*dt ^ are given by 

d Xl \dip = d0- ncrdt = 
d X2 \dip = d£-adt = 
d X3 \d V = 

<9 CT Jd?/>= (v / « 3 Te- Aik- A 2 )di = 
9 K JdV' = cr(K(K 2 + e)" 1/2 -A 1 )dt = ' 

-«9 e JdV> = dAi - A 2 dr/ + A 3 d£ = 
-%jdV> = dA 2 - A 3 d<9 = 
-9 J? Jd^ = dA 3 + A 2 d6' = 

The first three equations represent the horizontality restriction, the 4th en 5th equation represent 
the Euler-Lagrange optimization of the energy and the last three equations provide the Lagrange 
multipliers (recall (7.105 1) 

Ai = 



A 2 = -y/iy/l^ 

dz + A 3 erdi = dz + A 3 ds 



(A. 144) 



A 3 — — 7. 



by employing Noether's theorem and an invariance group of symmetries of ip (which is defined as 
a group acting on Q such that the induced action 77 on T*(Q) satisfies r) g {Z) = Z for all g € G) 
as will briefly explain next. In our case the action ij on T*(Q) is induced by the action of SE(2) 
acting on itself) of the minimization problem. 

Noether's theorem (which provides a conservation law on momentum) says that the momentum 
mapping m : Z — > SE(2)* given by 



(m(p) > = (CjV')(p), 



pez, 



is constant along the characteristic curves. Note that computation of the Lie-derivative of ij; along 
a characteristic curve gives 

= C^A = A\dijj + d{A\ip) = A\d%l), for all left-invariant vector fields A G C(SE(2)), 



which explains the last three equalities in (A. 143 1. 

The momentum mapping is invariant under the co-adjoint representation Ad* (this is the 
representation dual to the adjoint representation (6.52 1 ) 



m (Vg(p)) = (Ad -i)*m(p) 



(A. 145) 



which follows from the fact that r]*ip = ip and (?7 ff )*£ = (Ad g -i)*£. Consequently, the characteristic 
curves are contained in the co-adjoint orbits. It can be verified that the co-adjoint orbits of SE(2) 
are given by 



An + A 3 = c 2 e > 0, 



c>0, 



so we get the following preservation law that holds along the characteristic curves 



(z(s)) 2 + e-c 2 e = e(z( S )f, 



s > 0, 



where the normalized curvature z(s) — , K ^ ■■ 



(A. 146) 

satisfies \z\ < 1 and indeed this formula follows 



by integration of (7.1041, since 

(z(s)) 2 = e(z(s)) 2 + C,C € 



z = ez => zz = ezz 



As observed by Bryant Griffiths [9]p. 543-544 (with slightly different conventions) the last three 
equations of (A. 143) can be written 



dA = Xg^dg d(A • .g" 1 ) = d ((Ad 9 -i)* • A) = 
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where A = (—A3, A2, Ai) and where the matrix form of the Cartan-connection equals, 

/ t-os H - >lll H .r \ 

where both Lie-algebra and Lie-group are embedded in the group of invertible 3x3 matrices. 
Consequently, by Noether's theorem we have A = fx - g, for some constant fx = (—1x3,112,1x1), or 
more explicitly we have 

z = ni- li 3 x + Li 2 y 
z = —LI3 cos 9 + [12 sin 9 

\Je(l — z 2 ) = ll 3 sin 9 + li 2 cos 6, with fx 2 + Li 2 = c 2 e. 
Next we choose 



A*3 










c 2 e 


M2 










c 2 e 








1 



so that fx-hq 1 — ( v^c, 0, 0) and use left-invariance g = h^ 1 g , g = (x , y , e 10 ) then we get \ = fx-g 
(<Ve,0,0) • <7, i.e. 

i = ^tj, cy^cos = Cyfex — z, and — \J e(l — z 2 ) = — y^csin 9 = c^/ey 
and consequently we have 



x(s) = (y/ec) 1 z(s) 

9(s) = Z(i(s),e x ) = 9(0) + f k(t) dr. 





So by means of (7.1061 we get the solution g(s) = (x(s),y(s),9(s)) = h (x(s),y(s),0(s)), i.e. 



x(s) = ^ - ^x(s) - ^y(s) S( S ) = ^-cosh( v ^ S ) + 5sinh(7e S ) 

V(s) = + ^ e x(s) - ^y(s) with J y( s ) =y + l c J* V 1 ~ c 2 (x(r)) 2 e dr 

9(s) = 9(s) + arccos (-^g) I 9(s) = arccos Of sinh(^s) + cosh(^s)) , 

(A. 147) 



where c = \/l + z§. 



Now we have 6 unknown parameters fx%, fx 3 , z , z' Q , y(0), L, (note that \Xi is not unknown since 

, 2 J- i/2 — ^,2, j 1 ( z o) 



A* 2 + — c e an< i c = Y 1 + ( z o) 2 ) to ensure the given boundary conditions 

5 (0) = (x(0),y(0),e^)=g o := (*„, y , e ie °), 
g(L) = (x(L),y(L),e^) = g x = (x 1 ,y 1 ,e i ^) 

By means of left-invariance we can always make sure (by multiplying from the left with g^f 1 ) that 
g x = e, so 9x = 0, X\ = 0, y x = 0- 

In this case straightforward and intense computations yield 



y(o) 



Ml = z + - ^22/0, 

Li-2 = Cy/e sin (arccos (^fj), 

ll 3 = ~z' cos 9 + V~^in9oy/T^4, L={ % u* [ '"^'j^^™*'*"* ) a « > 1. w < o. 4 + *oV~< > 
_ I _ L (4) 2 7~7 2 



- 1 

sfi 




e 


3 

)« 


c = 1 


1 






V- 


-(4) 2 + (^0) 2e + C§ 


</i 








1 


log 


-M3- 


V- 









C - V e _ V + 6 V* ' ^: los V 7= ^ if c < 1, M3 < 0,4 + < 
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(A. 148) 



So all parameters are now expressed in the two unknown zq and z' which are determined by the 
two remaining boundary conditions: 



-»^ir(I)-Ml) 



Xi, 



c- < 



c-Je 



(A. 149) 



Now since SE(2) is a symmetric space [37] all points can be connected by a geodesic and we may 
expect that there indeed exist Zq and z' such that ( |A.149 1 holds. Consequently, the singularities 
(which cause extreme problems in the numerical shooting algorithm (7.951 of section 7.1l where 
z(s max ) — 1 occur always at s max > L (and if /i 3 ^ c^/e then s max > L). Next we explicitly verify 



that s r 



> L in 2 cases. 



In case c > 1 , 113 < and jz + ^/ezo > we have 

gV /IL _ -M3+-y/-(Zo) 2 +(zo) 2 e+M| 

and indeed — /X3 + \J — (z' ) 2 + (z ) 2 e + /Lt| < < (1 + c)-^ so L < s r , 
In case c < 1 , /i3 > and z' + < we have 



^yiL _ -M3-V-(z/)) 2 + (zo) 2 e+^| 

Zo+fWe 

^Vesmai. = -^+V / ( z o) 2 -( z o) 2£ + e 



ygO±£) 



and indeed we have e v/ " max > e^ L , since M3+ \/— (z ) 2 + ( z o) 2 e + A*| < Cyfe+ ^/e(l — c 2 ) + c 2 e 
^(1 + c). Equality is obtained if /i3 = Cy/e. 
See Figure [21] and see Figure [22] 
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Figure 22: Left: A phase plot of z(s) = K *^ +e and z(s) clearly indeed reveals that all paths 
will eventually end up at z — 1 where solutions brake down because of infinite curvature. Except 
for the cases where the initial condition is such that c = \J\ — e _1 (i(0)) 2 — (z(0)) 2 = 1 then 
solutions stay at c — 1 but reach the line z — 1 only for s — > oo. Right: These infinite curvature 
singularities always take place at s max > L, so this causes no problems in our exact analytic 



solutions (in contrast to the numerical shooting algorithm (7.951.) 



B Completion measures, Brownian bridges and geodesies 

on P Y = (SE{2),SE(2)/Y,tt,R) 

In this chapter we shall again work with the principal fiber bundle Py = (SE(2), SE(2)/Y, tt, R), 
ir(g) = [g] = gY, Y = {(0,y, 0) | y € K}, R g h — hg, equipped with Cartan connection ui : 
T(SE(2) -» T e (SE(2)) given by oj(X g ) = {it), X g )d x so that the horizontal part H of the tangent 
space T(SE(2) of SE{2) is (by definition) the kernel of the connection uj which is spanned by 

H = span{<9(3, d^}. 



Recall from section [673] that the tangent vectors along horizontal curves, recall Definition |5.5| are 
always in the horizontal part of the tangent space. 

Sofar we considered completion fields as collision probability of a forward direction process and 
a backward direction process. However, the direction process is a stochastic process for contour 
completion. In this section, however, we shall replace the (exact) resolvent Green's functions 
R^ 11 for contour completion by the Green's function of contour enhancement R^ 11,D22 . The 
corresponding completion distribution on SE(2) between e = (0,0, e l °) S SE(2) and h £ SE(2) 
is now given by 

g ^ R^ D ™{g)R^> D "{g- l h). (B.150) 

Next, we show it is related to the Brownian bridge measure by condition time integration, i.e. 
Laplace transform. 

The Brownian bridge measure on SE(2) is given by 

t 

A 

where A is a measurable set within SE(2) supported by the set of paths starting at time at e 
and ending up at h at time t > and where the density K.^ t : SE(2) — * M+ is given by 

VM= ' ^iH) (R151) 



GG 



where K^ 11 ^ 22 denotes the heat-kernel determined in subsection 4.1 The kernel Kg t represents 
the probability density that a random walker which started at e and ends up in h passes g. Note 
that it is a conditional measure, since it only considers paths with start at e and end at time t at 
h, which explains the denominator. We associate to such Brow nian br idge measure the following 



unconditional measure fC s ' t by removing the denominator in (|B.151|: 



K h s f c (g) = K^'^^K^'^ig-'h). 
Then we have the following result by Wittich, see [55] f° r more details and proof: 

Theorem B.13. Let M be a Riemannian manifold with distance function g?m : M — > R + with 
positive curvature. Let p £ M and choose r > 0, eo > such that R + £q < r(p), where r{p) > 
is some number such that B{p,r) := {x G M : \]d\{(p,x) < r} is strongly convex, that is for 
any two points within such a ball there is a unique minimizing geodesic whose interior is again 
contained within the ball. Let C C B(p,r) be an arbitrary closed subset. Write 

k(C) := max < 1, sup > , 

{ qeC,<ye/\ 2 (M)K(a) J 

where f\ q {M) denotes the set of anti- symmetric bilinear forms on the tangent space T q {M), which 
is isomorphic to the set of oriented 2D-planes in T q {M). Let Qq'' denote the Brownian Bridge 
measure on M supported by the set £l(p,q,t) of continuous paths starting at time at p and 
ending up at q G M at time t > and let Q t ' p ' q the associated unconditional measure, i.e. the 
conditional measure Qq'j and unconditional measure Q^' p ' 9 are given by 

OP<9(A\ _ i: m W P M (An{w(t)£B(q, v )}) 

w' M (Wt)eB(w)}) ' 
Q 7^(A) = lim W* f (A n {w(t) G B(q, r,)}) , 

where Wj, f denotes the well-known Wiener measure on M centered at p. Then for all q G C there 
is a unique geodesic jP^, parameterized by with constant velocity v(s) — dM , joining p and 
q. Furthermore, there is some 5 > such that for all e > with e < eo and for all q G C , 



%> q t {B{e,p,q,t))<2e 

®V?< M (B(e,p,q,t))<2e- c ~ 

where B(e,p,q,t) := {uj en(p,q,t) : sup d(u(s), ^'^(s)) > e} and R(n(C)) := VE^lSE^l cot ( 

se[o,t) \ 

Consequently, the Brownian Bridge measure tends, as 1 1 to the point measure supported by 
the geodesic ^p*^*, parameterized proportional to arc-length with constant velocity v(s) — dM (P' q > . 

Set M = SE(2) or rather M = Py, since we again restrict ourselves to horizontal curves. Then 
we set 

d S E(2)(9,h) = d SE(i) (e,g- l h) = inf { jj y/^( S ) + eds = J V(^'(*)) 2 + e||x'(t)|| 2 dt | 

7 is a smooth horizontal curve connecting g and go, 7(0) = a, 7(1) = g^h} , e = 

(B.152) 

where we recall ^ = ||x'(i)|j and we recall that for horizontal curves we have n(s) = 

Now we note that the ex plicit r elation between the completion measure fi e,h (A) induced by 
the completion distribution (B.150) and the unconditional Brownian-Bridge measure Q 7 ^' 6 '' 1 is 
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given by 

^ h {A) =jR^> D ™(g)R^> D "(g- 1 h)dv SBi2) (g) 



A 



a 2 / C(t ^ K^' u ™ (g))(a)C(t ^ K t u ^ (g- l h))(a) d^ SE{2) (g) 

A 

a 2 C (t^ / ^fK^^(g)K^ D ^(g) dsj dp SBW (g)) (a) 

j^jK^Hg)K^' D22 (g)dfi S E ( 2 ) (g) ] j d*) (a) 



a 2 £ It 

a 2 £(t^Q ^(>].Hn, 



Now we recall that the diffusion generator 9| + d$ is hypo-elliptic (it satisfies the Hormander 



condition, recall subsection 4.2.11 as a result the diffusion kernel satisfies similar estimates, that 



hold for Green's functions of elliptic operators: 



\Kt{g)\ < 



for some c > 0, [2] , [34] , [TT] , for more details on this particular case see Appendix [d] and conse- 
quently we have for A = i?(e, e, h, t) 

(j, e > h {A = B(e, e, h, t)) = a 2 C(t ^ Q^ e - q (A))(a) < a ce~ (d ^ ( '' h ^ +2(R - 5)r ' 2 f _ 2 e -at dtj 

= a 2 v ^g 1 ( % /(rfg J!;(2) ( e ,/ 1 ))^+2(K-a)e^) 
V( d s E (2)(e^)) 2 +2(i?-5)e 2 

withi?= (l/2)- v //«(C)d S£ ;(2)(e,/i) cot((l/2) y/K(C)d SE (2){e, /i)) and k(C) = sup K q (a)where K(cr) 

q£C 

is the horizontal curvature (i.e. the sectional curvature in the plane spanned by {dg,d^}. 

Consequently, we see that if the expected life time E(T) — ^ tends to the completion measure 
between two delta distributions 6 e and Sh tends to the point measure S^.h supported by the unique 

geodesic minimizing J \J k 2 + jjy^ ds (which we explicitly derived in section j 7" - 1 j ) connecting 
e and h. 

C "Snakes" in SE(2) based on completion fields of orienta- 
tion scores 

In this section we will formulate a variational problem, where the energy consists of two parts, 
with the goal of finding a sufficiently horizontal curve with given beginning and ending that fits the 
data C : SE(2) — > M + , where C denotes for example the completion distribution of an orientation 



score (3.181. Here the internal energy of such a curve is the elastica functional J K 2 (s)ds and the 
external energy of the curve, which takes care that the curve fits the data, is the total integral 
J e~ c ( J ( s ^ds, where s > denotes the arclength in ]R 2 of the projected curve x = PR27. So this 
is just a direct generalization of the variational methods for so-called snakes in image analysis 
where a curve in M 2 is supposed to fit the image data (x,y) <—> f(x,y) with given restrictions on 
the internal energy of the curve which usually consists of a length and curvature penalization. 
However, the "snakes" on SE(2) have in principle the advantage that they can deal with crossing 
contours. 

Next we will derive the corresponding Eulcr-Lagrange equation, but before we can continue 
we need a small result on differentiating a function SE(2) along a horizontal curve. 

Lemma C.14. A smooth curve in SE(2) given by s 1— > 'y(s) = £(s)ej(s) + i](s)e v (s) + 9(s)e$(s), 
with^(s) — (x(s),e ld ^) G SE{2), ej(s) = cos9(s)e x +sm9(s)e y , e n (s) = — sm9(s)e x +cos9(s)e y 



68 



and s > the arclength parameter of the projected curve x — Pr^T on the spatial plane, is a 
horizontal curve in SE(2) iff ^ = — £k. Moreover, for such curves we have -£ — KT] = 1. Now if 
we differentiate a smooth function on C : SE(2) — > M along a horizontal curve we get 

^Ch(s)) = d 6 C( 7 (s)) + K(8)doC(rf(a)) (C.153) 
Proof By straightforward differentiation we get 

7(s) - ^(s)e s {s)+r ] (s)e ri (s)+9(s)ee) = (f(«)-/s(«)»7(«))e c («) + (i7(s)+/c(a)C(a))e )J (a)+fl(«)e fl , 

which is horizontal iff 4 s = — £k in which case we have 1 = || ^| || = | ^| — k?7|. The sign follows 
from the restriction Z(x(s),e x ) = 9(s). Finally we note that by the chain-law we have 

£C( 7 («)) = (Q( 7 (s))de + C )7 (7*( S ))dr ? + C e (7*( S ))d^7(s)) 

Theorem C.15. Let g$,g\ € SE(2). Let L > be fixed. Consider a smooth positive function 
C : SE(2) -> M+ and dearie (7 : S , J B(2) -> M+ by C(g) = e~ c ^. Then a local maximum (or 
mode) of the following variational problem 

L 

arg min{£ l ( 7 ) = f(e~ c ^"^ + K 2 (s))||a;(s)|| ds | s i— > 'y(s) is a smooth horizontal curve with total length L 
o 

connecting 7(0) = go € SE(2) and j(L) — go <E SE(2), x— PR27} 

is a horizontal curve s i— > 7 *(s) = (cc*(s),e 4e ( s -*) which satisfies the following Euler-Lagrange 
equation: 

k 2 + 3k'(s) + C v (~f*(s)) - ^<%( 7 *(«)) = Xk(s), (C.154) 

where X is a Lagrange multiplier in H., due to </ie restriction to curves with length L. 
Proof Consider the following (normal) deviation on the horizontal curve 

7 *(s) = (x*( S ),0*( S )) ~ leS (s) = (x*(s)+e5(s)e n (s), 0* (s)+arctan 1 _ ),*> (C.155) 

with 5 smooth, compactly supported within (0,L), so 6(0) = 5(L) = S'(0) = S'(L) = and 

L L 

J K(s)5(s)ds = and J rc(s)(7( 7 *(s)) S(s)ds = (C.156) 


where s > represents the arc-length of t he curv e x* = P R 2 7*. Note that the total length of the 
deviated curve equals L + 0(e 2 ) since by (C.156 1 we have 

L 

n(s)5(s)ds = 

£ =o 

Moreover, this deviation yields a new horizontal curve "/ e $ with length L connecting go and g±. To 
this end we note that Z((l — eS(s)n(s))e^(s) + 8'(s)e v (s), e x ) — 9(s) + arctan ( jz^^^ p with 
k(s) — 6'(s). Note that by first order Taylor-expansion of the integrand around g(s) we get 



L 

d f \\±( x *(s) + e5(s) eri (s))\\ds 



lim 



l ™-J L \(c(r(s)) + e6(s)^\ +arctan{ 1 ^W }f I ) (1 - e5(s) K (s)) 
-C("/*(s))(k(s) + e5"(s) + e5(s) K 2 (s)) 2 (l - e5(a) K (s)) - k 2 (s) + 0(e 2 )} ds 

k(s)C(j*(s)) + 2k(s) + k 3 (s)^ ds 

+ 2k"(s) + k 3 (s) ) ds 




[ ac 




d 


SC 


I a?7 


7*(») 


ds 


as 


[ 




d 


ac 




7*(s) 


ds 


as 



7*0) 

f -'' -I" 'V' 

7*(s) 
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Then since 7* is the minimizer and since 6 satisfies the side conditions ( C. 156 1 the gradient of 
the Energy should be linear dependent on the gradients of the side conditions (which equal k and 
C(7*)k) and the Euler-Lagrange equations read: 



2k" (s) + k 3 (s) + 



dC 


_±\ 


dC 


\ 


drj 


ds 

7 *( s ) 




7*(s)J 



31 (s), 



(C.157) 



with qi(s) = Xik(s) + \ 2 C('j*(s))k(s), 



(C.158) 



where Ai and A2 are some Lagrange multipliers. We will show the final step A2 = 0, by systemat- 
ically checking for all horizontal curve pertubations. 

The same technique can be applied to horizontal curve deviations of the type 



(x*( S ),0*(s)) -> (x*(s) + e6(s)e £ (s), fl* (a) + arctan f^fl ) 



(C.159) 



where S is an arbitrary compactly supported smooth function within (0,L), but this is just a 
re-parametrization of the same curve and gives (C.153I. 



Finally, we consider horizontal curve deviations of the type 



(x(a), 0(a)) 1 ^ (x(a) + e<5(a)e„(a) 



k(s) 



(C.160) 



where S is an arbitrary compactly supported smooth function within (0, L) with J 5{s) da = 



which are again length preserving up to 0(e 2 ): 

L 



d e I ||-(xOO+e*(*K(*) 



and which yield 

K-\ S ) 



k(s) 



e ( (s))\\da 



, 8'{s)\' , XJ n S'(L) 8'(0) 
6(s)k(s)+ — (a)da = 0+- v 



e=0 



f/a 



7*W 



ds \ k(s) 



k(s) 



33(s) 



«(£) k(0) 



(C.161) 



wit h q.^j = A.^At(a) Now divide equation (C.153I by k and differentiate with respect to a > and 
by (C.157 1 we may substitute (Cg)' — C n — q\ + 2k'" + k 3 this yields the following equation 



d_ 
ds 



K-\ S ) 



dC 




dC 


_ ± | 






7*(«)} 




7 . w da \ 





q 1 (s)-2K"'(s) + K 3 (s) (C.162) 



from which we deduce <?i(a) = 93(a) and thereby A := Ai = A3 and A2 = from which the result 
follows. □ 

D Semigroups generated by subcoercive operators on SE(2) 



In this chapter we shall apply the general theory in |48j on weighted subcoercive operators on Lie 
groups, to our case of interest: The Forward Kolmogorov equation of the contour enhancement 



processes (2.151. From this general theory we will deduce that the closure of the generator of 



a contour enhancement process indeed generates a holomorphic semi-group with a smooth and 
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fast decaying Green's function on SE(2) (which we explicitly derived in section [4]). Moreover, 
we will derive Gaussian estimates for both exact and approximate kernel (derived in subsection 



4.2) and we will put our approach of approximation in a more general context yielding a contin- 



uous family of holomorphic semi-groups connecting the exact semigroup in subsection 4.1 and its 



Heisenberg" -approximation in subsection 4.2 We will show that the Gaussian estimates for the 



kernel a 1 R a (x, y, 6) are surprisingly sharp for the approximate resolvent case if a J. (i.e. infinite 



lifetime) which is exactly given by (4.401. This indicates that the Gaussian estimates of the exact 
kernel (nice for computations in the spatial domain) can be used as reasonable, somewhat rough, 
approximations of the exact convolution kernels if Dn << £>22- Furthermore, these Gaussian 
estimates can be used for taking regularized/Gaussian derivatives on orientation scores (similar 



to (4.31 I and (4.321 where D22 = -D33) for the (horizontal) case D33 = 0. 



Let G be a Lie group, with Lie Algebra T e {G) of dimension d with basis {A\, . . . , A d }- Let 

{A 1 ,...,A d ,}c{A u ...,A d }, d'<d, 

be an algebraic basis of the same Lie algebra, that is there exist an integer r (called the rank of 
the algebraic basis) such that 

0i :=span{A 1 ,...,A d '},g 2 ■= span{[fii,£|i]}, ... , g r = span{[fl r _i,0 r _i]} = T e (G). (D.163) 

Now for each element in A £ T e (G) there exists a minimum integer k £ N such that A £ 2k- We 
shall refer to this integer as the weight of A. In particular the weights of the basis elements Ai 
will be denoted by Wi for i = 1, . . . , d. Note that this particular convention of assigning weights, 
implies Wi = 1 iff i = 1, . . . , d! . We stress that this particular convention coincides with a special 
case of a reduced weighted algebraic basis {Ai, . . . A d '} in [58], where the filtration {0a}a>o should 
satisfy Q\ = {0} if A < 1, Q\ C fl„ for a ll A < fx and [g x , M ] C Sa+p for all A, \i > and Q r = 



for large r, need not be given by (D.163 1 and where the weights {wi}i = \ t .,^ d ' need not be equal to 
one, but should satisfy Ai ^ (J 0a for alH = 1, . . . , d! . 

\<Wi 

00 n 

Let J(d') denote the space of multi-indices associated to the algebraic basis, J{d!) — (J {1, ... , d'} K 

n=0 k=0 

and for all a = . . . , i n ) £ J(d') we associate the the Lie algebra element A a — Ai 1 ■ ■ ■ Ai n and 
the weighted length 

n 

Nl = 51 

k=l 

where n is the Euclidean length of a which we shall denote by |a| = n. If C : J(d') — » C is such 
that C(a) — if ||a|| > to, for some integer m £ N and if there exists an multi-index a £ J(d'), 
with 1 1 a 1 1 = to, such that C(a) 7^ 0, then C is called an m-th order form. To each m-th order form 
we associate the m-th order left invariant operator 

Ac := E C(a)A a 

ct€J(d'),ot=(ii,...,i n ) 

with A a = A il . ..Ai n = d7i(A fl ) . ..d1l{A in ). 

Definition D.16. Then C is said to be a G-weighted subcoercive form if £ 2N for all i = 

1, . . . , d' and there exist fi > 0, v £ R such that the Garding inequality holds 



Re{{4> 1 A c 4>)u(SE(2))} > H max \\A a \\u{SE(2)) 

\|a|<f 



Associated to group G and reduced weighted algebraic basis {Ai, . . . , A d '} one can construct 
a homogeneous Lie- Algebra Go, [44], with dilations (7t)t>o by means of 

[A,B] t = 1 r 1 (ht(A)MB)}) &7t([A,B] t ) = [ lt (A), lt (B)}, 

where j t (AA = t^A, for i = 1, . . . , d. 1 ' ' 
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Now (T e (G), [•, -]i=i) is the original Lie Algebra (T e (G), [•,-]) and (T e (G), [•, -] := [ V ]t|o) is a 
homogeneous Lie Algebra with dilations (74)00 which is uniquely determined by the filtration 
corresponding to the reduced algebraic basis. It can be shown that the reduced weighted algebraic 
basis {Ai, . . . , Ad'} is a reduced weighted algebraic basis for the Lie algebra (T e (G), [•, •]<) for all 
t > 0. The group simply connected group Gt is generated by the Lie algebra (T e (G), [■, -]t) via the 
exponential mapping. The left-invariant vector fields A\ on Gt are given by 



(A\<t>){g) 



(d1Z(Ai)ct>)(g) = — <t>(gexp t {3Ai)) 
as 



Now the Lie algebra (T e (G), [•, -]t) can be equipped with the following modulus 
\g\' t = d' t (g,e)=mf{6>0 : 3 ieCt(S) : 7 (Q) = e, 7 (l) = g] , 



(D.166) 



where C t {5) equals the space of all absolutely continuous curves 7 with tangent vectors in the 
plane spanned by {A\ , . . . , A l d , } such that 



7(a) 



7i(«) ^| 7(>) , with = |<<L4*| tW ,t(*)>1 < foralU = l,...,d'and S >0. 



If we return to our special choice of filtration (D.163) all weights of the algebraic basis elements 
are equal 1 in which case (T e (G), [•, -]o) is a nilpotent Lie algebra of the same rank r, [35] Lemma 
3.10, p. 106. Using (T e (G), [•, -]o) as a "local approximation" of (T e (G), [-, ■]) the authors in [15] 
obtained their general result [15]Thm 1.1. p. 93. Next we give a brief summary for the special case 
of our interest (H,G,U) = (L2 (G) , G, 1Z) , G — SE{2), with 1Z the right regular representation 
whose derivative is the isomorphism between T e {G) and the Lie algebra of left invariant vector 



fields C(G), recall (2.7) 



Theorem D.17. Let C be a G-weighted subcoercive form defined on a Lie group G, with Haarmea- 
sure hg- Then the closure of —Ac generates a holomorphic semigroup s t— > S s on L,2(G) which 
has a fast decreasing kernel in K s € Li(S'£ , (2)) D C°°{G) such that 



A a S s U = / (A a Ks)^- 1 g)U(h)dfi G (h), for all a e J(d'), and all U e L 2 (G) 



(D.167) 



and for all a there exist b, c > such that 



\A a K s (g)\ <cs 



'" "(l|a|| + E -6' — 

=1 e 



(D.168) 



for all g G G and all s > 0, where \g\± is given by (D.166) with t = 1 



D.l Application to Forward Kolmogorov Equation of Contour Enhance- 
ment Process 

Now to apply this result to the case of the Forward Kolmogorov equation for the contour enhance- 
ment process we set 

G = SE(2),T e {SE(2)) = span^, A 2 , A 3 } = speai{de,d x ,d y }, 
C(SE(2)) = spim{Ai,A 2 ,A 3 } = span{9 e , d ( , d v }. 

As algebraic basis we use {Ai,^} since the corresponding left-invariant vector fields {^1,^2} 
span the horizontal part of the tangentspaces with respect to the Cartan-Maurer form ( 6.54| , 



which is the only natural choice if it comes to the restriction of curves in SE{2) which arise as 
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"lifts" from their projections on the spatial plane. The rank of this algebraic basis is equal to 2 
and we have 

0i = span{<9 e ,<%} and g 2 = [fli,0i] = span{<9 e , d e , d^} = £(SE(2)), 



so Wi — W2 — 1 and d v S 02 => w 3 = 2. Now C(a) — —a\ — a§ is a SE{2) weighted subcoercive 
form, since m — 2 and U)i = W2 = 1 and moreover the Garding inequality holds with /i = 1, v = {) 
since by partial integration 

-(0, (<9^ + <9f )0)l 2 (S£(2)) > max{||9e0||L 3 (S£;(2)). II 5 ^IIl 2 (S£;(2))} 

for all test functions <f> S 2?(f2 
Consequently, by Theorem 



D.17 



where fl e is some open environment around the unity in SE{2). 
the closure of d 2 + <9| generates a holomorphic semigroup t ^ St 



on L 2 (S'i;(2)) with a fast decreasing kernel K t e h 2 (SE(2)) n Li(5S(2)) such that ( |D.167| | holds 
Moreover, the kernel satisfies the following estimate: 

\K t { 9 )\<ct-*e-^, 

with locally (|<?|') = £ 2 + 9 2 + \rj\, since (0, £, n) are the coordinates of the second kind in SE(2). 
The kernel satisfies (d t + d 2 + d^)K t = 5$ ® <5f as distributions. 

This can be generalized to D\\ > 0, D22 > and C(a) = —Duct 2 — D 22 o!\ yeilding 



K t (x,y,0) < 



1 



47Ti 2 L> 11J D 2 2 



V /D H- D 22 



'}* 



for all x, y € R and all 9 £ [0,27r), where we note that the constant b in the equality (D.168I is 
independent on Dn,D 22 so b — j. This estimate also coincides with the estimate by Citti and 
Sarti [TTJ Thm 5.1. 

Next we estimate will investigate the sharpness of the Gaussian estimates. Now since the 
moduli I • \ t are locally equivalent [48] and since we have a simple exact formula for the resolvent 



Heisenberg approximation kernel (4.40) in the spatial domain we choose to study the sharpness of 



the estimate of the kernels on the group Go which will turn out to be isomorphic to the Heisenberg 
group H3. Here we shall devide the analysis in two steps. First we shall show that the approach 



in subsection 4.2 is a special case of the above homogcnization (D.165) of the Lie algebra, yielding 



a continuum of semigroups between the exact case studied in subsection |4.1| and the Heisenberg 
approximation case studied in subsection |4.2| Then we shall consider the logarithmic weighted 



modulus which is locally equivalent to the modulus (D.166 ) and derive surprisingly sharp Gaussian 



estimates both from above and below of the resolvent Heisenberg kernel (4.401 



Following the general scheme we define the dilation on the algebra by jt '■ T e (SE(2)) — » 
T e (SE{2)) by 7t(c 1 Ai + c 2 A 2 + c 3 A 3 ) = tc l Ax + tc 2 A 2 + t 2 c 3 A 3 . Furthermore we define the 
corresponding dilation on the group by j t (x,y,e te ) — (f, J^e 1 *). Now (T e (SE(2)),[-,-] t ) with 



[•,•]* given by (D.165 1 is a Lie-algebra with corresponding simply connected group (SE(2)) t = 
exp t (T e (SE(2))). Note that the dilation on the Lie-algebra coincides with the pushforward of the 
dilation on the group j t = (7)* and the left invariant vector fields on (SE(2)) t are given by 



A\\ g ={% °L g oj t )*Ai, 
for all t € (0, 1] and a brief computation yields 



A- 1 



{% L g)*{lt)*Ai(f) ■■ 

t 1 J\a 



(% 1 Lg)*7t(Ai 



= f 



'{it l 0Lg)*(Ai 



= f 



'(ir 1 ). 



l % -^(0°7t) 



for all smooth complex- valued functions <fi defined on a small open environment around g € SE{2). 
So we see that for all g = (x, y, e l6 ) 6 SE(2) we have 



A\\ 



Hide) = d e 



cos(gt) q _|_ gin m . 
f sin(0£) q cos(dt) ^ 



cos(6t)d x 



i(6t) 



On 



-tsm(0t)d x + cos(6t)d y 
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and indeed A£\„ = Ai\„ = A; and 



[A\,A\] = 0, [A\,A\] = t 2 A\, [A{,A 3 ] = t 2 A\, 

and by taking the limit t j we see that the homogeneous contraction (SE(2))q — lim(SE(2)) t 

is isomorphic to H 3 and the space of corresponding left-invariant vector fields equals C(H 3 ) = 
span{<9#, d x + 8d y , d y } which is indeed a nilpotent Lie group of rank 2. This nilpotent Lie group 
isomorphic to H 3 is a subgroup of the five dimensional group H§ of Hciscnbcrg type that arises 
by approximating cos 8 f=s 1 and sin 8^8 whose left-invariant vector fields are given by 

A 1 = Al = d e , M = Al = d x + 6d y , A 4 = d y , 
A 3 = -8d x + d v , A 5 = d x . 

This Lie-algebra C(H§) = span{Ai, A\, A3, A4, A§) is isomorphic to the matrix-algebra 



i=i 

whose exponent is given by 
exp(tB) = 1 + tB 



( a 1 a 4 a 5 \ 

a 2 a 3 



\ / 



B 



i=l 



f 2 

-rB 2 = 



1 


ta 1 


ta 4 + \t 2 a}a 2 


to 5 + \t 2 a x a 3 





1 


ta 2 


ta 3 








1 














1 



This isomorphism enables us to quickly relate the coordinates of first kind to the coordinates of 
the second kind in H3 = (SE(2))q without explicitly using the CBH-formula : 

(a;, y, 8) = exp (a 3 A 3 ) exp (a 2 A 2 ) exp(a 1 Ai) = exp (/3Mi + (5 2 A 2 + f3 3 A 3 ) <^ 



/3 1 = a 1 = 8, (3 3 + i/3 1 /? 2 = a 3 = y, f3 2 



so we see that the coordinates of the first kind on (SE(2))o read 
f] 1 = 9, f3 2 = x and f3 3 = y - ]-x6 



and as a result the weighted modulus on (SE(2))q associated to the filtration (D.163I is given by 



\g\ = \j8 2 + x 2 + \y--x8\, 

now by @H] Prop. 6.1 there exists a c > 1 and an e S (0,1] such that for all a e T e (SE(2)) 
with ||a|| < e such that c _1 |a| t < |exp t (a)|J < c\a\t, where the weighted modulus is given by 

\a\ t = \J2Pl A l\ = V(Pl) 2/wi + WF™ 2 + \Pt\ 2/w2 - As a result we have for t = 0, $L = (3 k , 
k = 1,2,3 that 



\(x,y,8)\'>l(x 2 + 8 2 + \y-lx8\)^ 
kf^{x,y^)<^e- i ^^ < 



1 P -JrA x + e +\y-k^\) 

_„2 e <= * s , 



so that by integration over traveling time s > we find 

00 

]imR? lx > D ™(x,y,6) = [ kf luD22 (x, y, 8) ds < 1 



D22 jsL + 
D22 T Oil 



\/-Dll-D22 
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Now if we consider the exact solution 

< J K^' D ™(x,y,e) As 



< 1 75 

— 1D11C22 m 2 o2 l«-4'" 
-t>22 + Oll + v/Dn D 22 



where we note that for all a, 6 > one has a + 6 > V« 2 + fr 2 > ^75 ( a + then we sec that 

c = \/2 w 1.19 > 1 indeed yields a Gaussian upper bound for the exact Heisenberg kernel for 
a I 0, whereas c = 0.5 yields a Gaussian lower-bound for the same kernel. 
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